In Learning from Demonstration (LfD), a human (or another controller) provides example executions of a task. Behavior cloning directly regresses a policy from states to actions, but it does not explicitly model why an action is good, and it often suffers from covariate shift when the robot leaves the demonstrated state distribution.

Inverse reinforcement learning instead assumes demonstrations come from (near-)optimal behavior under some unknown reward function \( R^\star \). The expert policy \( \pi_E \) is assumed to satisfy

\[ \pi_E \approx \arg\max_{\pi} J(\pi; R^\star), \]

where \( J(\pi; R) \) is the expected discounted return of policy \( \pi \) under reward \( R \). IRL seeks to reconstruct a reward \( R_\theta \) that makes the expert appear optimal, and then a standard RL algorithm (to be studied systematically in Chapter 12) can be used to compute a robot policy for new environments.

For robotics, this is compelling because a reward function can generalize to new initial conditions, obstacles, or robot morphologies, while raw action demonstrations often cannot.

flowchart TD
  D["Expert demonstrations (state-action sequences)"] --> F["Estimate feature statistics of expert behavior"]
  F --> IRL["Solve inverse RL optimization to find reward parameters"]
  IRL --> RWD["Recovered reward function R_theta"]
  RWD --> RL["Forward RL planning / control with R_theta"]
  RL --> ROBOT["Robot executes derived policy"]
        

We model the robot and its environment as a Markov Decision Process (MDP):

\[ \mathcal{M} = (\mathcal{S}, \mathcal{A}, P, \gamma, \rho_0), \]

where \( \mathcal{S} \) is the (possibly high-dimensional) state space (e.g., joint configurations and object poses), \( \mathcal{A} \) is the action set (e.g., end-effector twists, joint torques, or motion primitives), \( P(s' \mid s, a) \) is the transition kernel, \( \gamma \in [0,1) \) is the discount factor, and \( \rho_0 \) is the distribution of the initial state.

For a stationary policy \( \pi(a \mid s) \), the discounted return under reward \( R : \mathcal{S} \times \mathcal{A} \to \mathbb{R} \) is

\[ J(\pi; R) \;=\; \mathbb{E}\!\left[ \sum_{t=0}^{\infty} \gamma^t R(s_t, a_t) \right]. \]

A key quantity in IRL is the discounted occupancy measure of a policy:

\[ d_\pi(s,a) \;=\; (1 - \gamma) \sum_{t=0}^{\infty} \gamma^t \Pr_\pi(s_t = s, a_t = a). \]

Proposition. For any bounded reward function \( R \), we have

\[ J(\pi; R) \;=\; \frac{1}{1 - \gamma} \sum_{s \in \mathcal{S}} \sum_{a \in \mathcal{A}} d_\pi(s,a) \, R(s,a). \]

Proof. Starting from the definition,

\[ \begin{aligned} J(\pi; R) &= \mathbb{E}\!\left[ \sum_{t=0}^{\infty} \gamma^t R(s_t, a_t) \right] \\ &= \sum_{t=0}^{\infty} \gamma^t \sum_{s,a} \Pr_\pi(s_t = s, a_t = a)\, R(s,a) \\ &= \sum_{s,a} R(s,a) \sum_{t=0}^{\infty} \gamma^t \Pr_\pi(s_t = s, a_t = a) \\ &= \frac{1}{1 - \gamma} \sum_{s,a} d_\pi(s,a) R(s,a), \end{aligned} \]

where the last line uses the definition of \( d_\pi(s,a) \). This representation is central because IRL often works in terms of occupancy measures rather than policies directly.

Let \( \pi_E \) be the expert policy that generated the demonstrations and \( d_E(s,a) = d_{\pi_E}(s,a) \) its occupancy measure. The IRL problem assumes \( P, \gamma, \rho_0 \) are known (or estimated) and that the reward has a parametric form:

\[ R_\theta(s,a) = \boldsymbol{\theta}^\top \boldsymbol{\phi}(s,a), \quad \boldsymbol{\theta} \in \mathbb{R}^k, \]

where \( \boldsymbol{\phi}(s,a) \in \mathbb{R}^k \) is a feature vector (e.g., distances to obstacles, joint limit penalties, manipulability indices). From Section 2,

\[ J(\pi; R_\theta) = \frac{1}{1 - \gamma} \sum_{s,a} d_\pi(s,a)\,\boldsymbol{\theta}^\top \boldsymbol{\phi}(s,a) = \frac{1}{1 - \gamma} \boldsymbol{\theta}^\top \boldsymbol{\mu}(\pi), \]

where the feature expectation of policy \( \pi \) is

\[ \boldsymbol{\mu}(\pi) = \sum_{s,a} d_\pi(s,a)\,\boldsymbol{\phi}(s,a) = \mathbb{E}_\pi \left[ \sum_{t=0}^{\infty} \gamma^t \boldsymbol{\phi}(s_t,a_t) \right]. \]

Given expert trajectories \( \{\tau^{(i)}\}_{i=1}^N \), we estimate \( \boldsymbol{\mu}_E \approx \boldsymbol{\mu}(\pi_E) \) via

\[ \hat{\boldsymbol{\mu}}_E = \frac{1}{N} \sum_{i=1}^N \sum_{t=0}^{T_i - 1} \gamma^t \boldsymbol{\phi}\!\left( s_t^{(i)}, a_t^{(i)} \right). \]

A generic IRL objective is to find a reward such that the expert performs at least as well as any other policy:

\[ \boldsymbol{\theta}^\star \in \arg\max_{\boldsymbol{\theta}} \Big( J(\pi_E; R_{\boldsymbol{\theta}}) - \max_{\pi} J(\pi; R_{\boldsymbol{\theta}}) \Big). \]

This is intractable in general, but it motivates practical formulations such as max-margin IRL and maximum entropy IRL, which we discuss next.

A classic IRL approach for robotics is apprenticeship learning via feature expectation matching. We maintain a set of candidate policies \( \{\pi_1,\dots,\pi_K\} \) and their feature expectations \( \boldsymbol{\mu}(\pi_j) \). The reward parameters are chosen to maximize the margin between the expert and all current candidates:

\[ \begin{aligned} &\max_{\boldsymbol{\theta}, \delta} \;\; \delta \\ \text{s.t.}\quad &\boldsymbol{\theta}^\top \bigl( \hat{\boldsymbol{\mu}}_E - \boldsymbol{\mu}(\pi_j) \bigr) \ge \delta, \quad j = 1,\dots,K, \\ &\|\boldsymbol{\theta}\|_2 \le 1. \end{aligned} \]

Claim. The optimization above is a convex program (a second-order cone program).

Proof sketch. The constraints \( \boldsymbol{\theta}^\top(\hat{\boldsymbol{\mu}}_E - \boldsymbol{\mu}(\pi_j)) \ge \delta \) are linear in the decision variables \( (\boldsymbol{\theta},\delta) \), and the norm constraint \( \|\boldsymbol{\theta}\|_2 \le 1 \) is convex. The feasible set is therefore convex, and the objective \( \delta \) is linear. Maximizing a linear function over a convex set is a convex optimization problem.

Once \( \boldsymbol{\theta} \) is obtained, we solve a forward RL problem under reward \( R_{\boldsymbol{\theta}} \) to compute \( \pi_{K+1} \). If its feature expectations are already close to \( \hat{\boldsymbol{\mu}}_E \), the procedure stops; otherwise, we add \( \pi_{K+1} \) to the pool and resolve the max-margin problem. This yields a sequence of policies that converges (under suitable assumptions) to one whose feature expectations are close to those of the expert.

Max-margin IRL returns a reward that separates expert and non-expert feature expectations, but it does not specify a full distribution over trajectories. Maximum entropy IRL (MaxEnt IRL) instead chooses the most uninformative distribution consistent with the expert feature expectations.

Let \( \mathbf{f}(\tau) \) denote the discounted feature sum along a trajectory:

\[ \mathbf{f}(\tau) = \sum_{t=0}^{T(\tau)-1} \gamma^t \boldsymbol{\phi}(s_t, a_t). \]

MaxEnt IRL posits a Gibbs distribution over trajectories:

\[ p_{\boldsymbol{\theta}}(\tau) = \frac{1}{Z(\boldsymbol{\theta})} \exp\bigl(\boldsymbol{\theta}^\top \mathbf{f}(\tau)\bigr), \]

where \( Z(\boldsymbol{\theta}) \) is the partition function. Given demonstrations \( \{\tau^{(i)}\}_{i=1}^N \), the log-likelihood is

\[ \mathcal{L}(\boldsymbol{\theta}) = \sum_{i=1}^{N} \log p_{\boldsymbol{\theta}}(\tau^{(i)}) = \sum_{i=1}^{N} \boldsymbol{\theta}^\top \mathbf{f}(\tau^{(i)}) - N \log Z(\boldsymbol{\theta}). \]

Proposition. The gradient of the log-likelihood is the difference between expert and model feature expectations:

\[ \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) = \hat{\boldsymbol{\mu}}_E - \boldsymbol{\mu}_{\boldsymbol{\theta}}, \]

where

\[ \boldsymbol{\mu}_{\boldsymbol{\theta}} = \mathbb{E}_{p_{\boldsymbol{\theta}}} \bigl[\mathbf{f}(\tau)\bigr]. \]

Proof. We compute

\[ \begin{aligned} \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) &= \sum_{i=1}^{N} \mathbf{f}(\tau^{(i)}) - N \frac{1}{Z(\boldsymbol{\theta})} \nabla_{\boldsymbol{\theta}} Z(\boldsymbol{\theta}) \\ &= \sum_{i=1}^{N} \mathbf{f}(\tau^{(i)}) - N \sum_{\tau} \mathbf{f}(\tau) \frac{ \exp\bigl(\boldsymbol{\theta}^\top \mathbf{f}(\tau)\bigr) }{ Z(\boldsymbol{\theta}) } \\ &= N \hat{\boldsymbol{\mu}}_E - N \mathbb{E}_{p_{\boldsymbol{\theta}}} \bigl[\mathbf{f}(\tau)\bigr] = N\bigl( \hat{\boldsymbol{\mu}}_E - \boldsymbol{\mu}_{\boldsymbol{\theta}} \bigr). \end{aligned} \]

Dividing by \( N \) gives the stated expression. Gradient ascent iteratively adjusts \( \boldsymbol{\theta} \) to align model and expert feature expectations.

In practice, we do not enumerate trajectories. Instead, MaxEnt IRL can be implemented via soft value iteration, where the soft value and soft Q-functions satisfy

\[ Q_{\boldsymbol{\theta}}(s,a) = R_{\boldsymbol{\theta}}(s,a) + \gamma \sum_{s'} P(s' \mid s,a) V_{\boldsymbol{\theta}}(s'), \]

\[ V_{\boldsymbol{\theta}}(s) = \log \sum_{a \in \mathcal{A}} \exp\bigl(Q_{\boldsymbol{\theta}}(s,a)\bigr), \]

and the induced stochastic policy is

\[ \pi_{\boldsymbol{\theta}}(a \mid s) = \exp\bigl(Q_{\boldsymbol{\theta}}(s,a) - V_{\boldsymbol{\theta}}(s)\bigr). \]

Most IRL algorithms used in robotics (for manipulators or mobile robots) follow an iterative loop that alternates between forward planning and inverse reward updating:

flowchart TD
  S["Start with expert demonstrations"] --> MU["Estimate expert feature expectation mu_E"]
  MU --> THETA0["Initialize reward parameters theta_0"]
  THETA0 --> LOOP["Repeat: forward RL + reward update"]
  LOOP --> FWD["Forward RL: solve MDP with current R_theta to get policy pi_k"]
  FWD --> MU_K["Compute mu(pi_k) from induced trajectories"]
  MU_K --> UPD["Update theta_{k+1} using mu_E - mu(pi_k)"]
  UPD --> CHECK["Check convergence in feature space"]
  CHECK -->|not converged| LOOP
  CHECK -->|converged| OUT["Return R_theta and policy for deployment"]
        

In a robotics system, the forward RL step is often implemented using continuous-control solvers (e.g., trajectory optimization, model-predictive control, or policy gradient methods) while IRL provides a task-specific cost shaping signal that can generalize to new contexts.

We now sketch a compact Python implementation of MaxEnt IRL on a discrete 2D configuration grid (e.g., a coarse approximation of an end-effector workspace). For real robotic systems, one would typically use numpy, scipy, and robotics simulators such as pybullet or mujoco to generate demonstrations, but the core IRL logic is similar.

import numpy as np

# Small gridworld as proxy for a low-res C-space
side = 3
n_states = side * side
n_actions = 4  # 0: up, 1: right, 2: down, 3: left
gamma = 0.9

def state_to_rc(s):
    return divmod(s, side)

def rc_to_state(r, c):
    return r * side + c

def make_transition():
    P = np.zeros((n_states, n_actions, n_states))
    for s in range(n_states):
        r, c = state_to_rc(s)
        for a in range(n_actions):
            nr, nc = r, c
            if a == 0:
                nr = max(r - 1, 0)
            elif a == 1:
                nc = min(c + 1, side - 1)
            elif a == 2:
                nr = min(r + 1, side - 1)
            else:
                nc = max(c - 1, 0)
            s_next = rc_to_state(nr, nc)
            P[s, a, s_next] = 1.0
    return P

P = make_transition()

# Feature map: one-hot over states (k = n_states)
k = n_states
phi = np.eye(n_states)  # phi[s, :] is feature vector for state s

def empirical_feature_counts(demos):
    """demos: list of trajectories, each as list of states"""
    mu_E = np.zeros(k)
    for traj in demos:
        for t, s in enumerate(traj):
            mu_E += (gamma ** t) * phi[s]
    mu_E /= len(demos)
    return mu_E

def soft_value_iteration(theta, n_iters=50):
    # reward only depends on state here: R(s) = theta^T phi(s) = theta[s]
    R = theta
    V = np.zeros(n_states)
    for _ in range(n_iters):
        Q = np.zeros((n_states, n_actions))
        for s in range(n_states):
            for a in range(n_actions):
                Q[s, a] = R[s] + gamma * np.dot(P[s, a], V)
        # soft value function
        V = np.log(np.exp(Q).sum(axis=1) + 1e-8)
    # derive policy
    Q = np.zeros((n_states, n_actions))
    for s in range(n_states):
        for a in range(n_actions):
            Q[s, a] = R[s] + gamma * np.dot(P[s, a], V)
    pi = np.exp(Q - V[:, None])
    pi /= pi.sum(axis=1, keepdims=True)
    return V, pi

def state_visitation(pi, start_dist, horizon=20):
    d_s = np.zeros(n_states)
    d = start_dist.copy()
    for t in range(horizon):
        d_s += (gamma ** t) * d
        d_next = np.zeros(n_states)
        for s in range(n_states):
            for a in range(n_actions):
                for s2 in range(n_states):
                    d_next[s2] += d[s] * pi[s, a] * P[s, a, s2]
        d = d_next
    d_s *= (1.0 - gamma)
    return d_s

def expected_features(theta, start_dist):
    _, pi = soft_value_iteration(theta)
    d_s = state_visitation(pi, start_dist)
    return d_s @ phi  # (n_states,) @ (n_states, k) = (k,)

# Example expert demonstrations (hand-designed "right then down" strategy)
demos = [
    [0, 1, 2, 5, 8],   # trajectory 1
    [1, 2, 5, 8],      # trajectory 2
]

mu_E = empirical_feature_counts(demos)

theta = np.zeros(k)
start_dist = np.zeros(n_states)
start_dist[0] = 1.0  # fixed start state

alpha = 0.5
for itr in range(50):
    mu_theta = expected_features(theta, start_dist)
    grad = mu_E - mu_theta
    theta += alpha * grad  # gradient ascent on MaxEnt IRL objective

print("Learned reward (per state):")
print(theta.reshape((side, side)))
      

In a manipulator setting, one could choose features that depend on joint configurations, distances between end-effector and target, and collision penalties obtained from a motion planning library (e.g., OMPL via Python bindings). The IRL loop above remains unchanged; only the feature extraction and transition model become more sophisticated.

In C++, IRL for robotics is often integrated with ROS, OMPL, and MoveIt. Below is a minimal snippet computing discounted feature expectations for demonstrations using Eigen for linear algebra. In a real system, phi(s,a) would call into robot kinematics libraries such as orocos_kdl or MoveIt to compute task-relevant features.

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

using Eigen::VectorXd;
using Eigen::MatrixXd;

struct Transition {
    int s;
    int a;
};

struct Demonstration {
    std::vector<Transition> traj;
};

VectorXd phi(int s, int a, int k);

// Empirical discounted feature expectation from demonstrations
VectorXd empiricalFeatureExpectation(const std::vector<Demonstration>& demos,
                                     double gamma,
                                     int k) {
    VectorXd mu_E = VectorXd::Zero(k);
    for (const auto& demo : demos) {
        double discount = 1.0;
        for (std::size_t t = 0; t < demo.traj.size(); ++t) {
            int s = demo.traj[t].s;
            int a = demo.traj[t].a;
            mu_E += discount * phi(s, a, k);
            discount *= gamma;
        }
    }
    mu_E /= static_cast<double>(demos.size());
    return mu_E;
}

// Simple gradient update for theta in linear reward R_theta(s,a) = theta^T phi(s,a)
void gradientUpdate(VectorXd& theta,
                    const VectorXd& mu_E,
                    const VectorXd& mu_theta,
                    double alpha) {
    VectorXd grad = mu_E - mu_theta;
    theta += alpha * grad;
}
      

Forward RL in C++ can be implemented using dynamic programming on discretized configuration spaces, or by delegating to trajectory optimizers (e.g., CHOMP, STOMP) and interpreting their costs as approximations to \( J(\pi; R_\theta) \).

Java-based robotic systems (e.g., via ROSJava or custom industrial APIs) can incorporate IRL by maintaining feature expectations and updating reward parameters online. The following snippet shows a basic structure:

import java.util.List;

public class LinearIrl {

    // phi(s,a) is assumed known and returns double[k]
    public static double[] phi(int s, int a, int k) {
        double[] f = new double[k];
        // TODO: fill with task-specific features
        return f;
    }

    public static double[] empiricalFeatureExpectation(List<int[][]> demos,
                                                       double gamma, int k) {
        double[] muE = new double[k];
        for (int[][] traj : demos) {
            double discount = 1.0;
            for (int t = 0; t < traj.length; t++) {
                int s = traj[t][0];
                int a = traj[t][1];
                double[] f = phi(s, a, k);
                for (int i = 0; i < k; i++) {
                    muE[i] += discount * f[i];
                }
                discount *= gamma;
            }
        }
        int n = demos.size();
        for (int i = 0; i < k; i++) {
            muE[i] /= (double) n;
        }
        return muE;
    }

    public static void gradientUpdate(double[] theta,
                                      double[] muE,
                                      double[] muTheta,
                                      double alpha) {
        for (int i = 0; i < theta.length; i++) {
            double grad = muE[i] - muTheta[i];
            theta[i] += alpha * grad;
        }
    }
}
      

A forward planner (e.g., a lattice-based planner for mobile robots) can be treated as a black-box policy generator under \( R_\theta \), from which Java code collects trajectory statistics \( \boldsymbol{\mu}(\pi_\theta) \).

MATLAB with the Reinforcement Learning Toolbox and Robotics System Toolbox is widely used for robotic experiments. Below is a simplified script implementing soft value iteration under a linear reward. The resulting policy can be wrapped in a Simulink block to close the loop with a simulated robot plant.

function [V, pi] = soft_value_iteration(theta, P, gamma, nIter)
% theta: k-by-1 reward parameters
% P: nS-by-nA-by-nS transition probabilities
% Here we assume state features are one-hot, so R(s) = theta(s).

nS = size(P, 1);
nA = size(P, 2);

R = theta(:);          % nS-by-1
V = zeros(nS, 1);

for it = 1:nIter
    Q = zeros(nS, nA);
    for s = 1:nS
        for a = 1:nA
            Q(s,a) = R(s) + gamma * squeeze(P(s,a,:))' * V;
        end
    end
    V = log(sum(exp(Q), 2) + 1e-8);
end

Q = zeros(nS, nA);
for s = 1:nS
    for a = 1:nA
        Q(s,a) = R(s) + gamma * squeeze(P(s,a,:))' * V;
    end
end

pi = exp(Q - V);
pi = pi ./ sum(pi, 2);
end
      

In Simulink, one can implement this as a MATLAB Function block that receives the current state and outputs a distribution over discrete actions, or is combined with a low-level controller that maps actions to joint torques.

Wolfram Mathematica is convenient for prototyping IRL with symbolic and numerical tools. The following snippet performs a few gradient ascent iterations on the MaxEnt IRL objective given empirical feature expectations and a function that returns model feature expectations.

(* Empirical feature expectation from demonstrations *)
gamma = 0.9;
demos = {
  {0, 1, 2, 5, 8},
  {1, 2, 5, 8}
};

phi[s_] := UnitVector[9, s + 1]; (* one-hot over 9 states *)

muE =
  Mean@Table[
    Sum[gamma^t phi[traj[[t + 1]]], {t, 0, Length[traj] - 1}],
    {traj, demos}
  ];

(* Placeholder for model feature expectation given theta *)
muTheta[theta_] := (* compute via soft value iteration and visitation *)
  Module[{(* ... implementation omitted ... *)},
    ConstantArray[0., Length[theta]]
  ];

alpha = 0.5;
theta = ConstantArray[0., 9];

Do[
  grad = muE - muTheta[theta];
  theta = theta + alpha grad;
, {itr, 1, 20}];

thetaMatrix = Partition[theta, 3]
      

In more advanced applications, Mathematica's NDSolve and control systems functions can be used to simulate robot dynamics while IRL updates the reward shaping terms.

Problem 1 (Occupancy Measure Representation). Let \( d_\pi(s,a) \) be the discounted occupancy measure as in Section 2. Prove that for any bounded reward \( R \), we have \( J(\pi; R) = \frac{1}{1-\gamma} \sum_{s,a} d_\pi(s,a) R(s,a) \).

Solution. By definition,

\[ J(\pi; R) = \mathbb{E} \left[ \sum_{t=0}^{\infty} \gamma^t R(s_t,a_t) \right] = \sum_{t=0}^{\infty} \gamma^t \sum_{s,a} \Pr_\pi(s_t = s,a_t = a) R(s,a). \]

Rearranging the sums gives

\[ J(\pi; R) = \sum_{s,a} R(s,a) \sum_{t=0}^{\infty} \gamma^t \Pr_\pi(s_t = s,a_t = a). \]

Multiplying and dividing by \( 1-\gamma \), and using the definition of \( d_\pi(s,a) \), yields the desired result.

Problem 2 (Linear Reward and Feature Expectations). Suppose \( R_\theta(s,a) = \boldsymbol{\theta}^\top \boldsymbol{\phi}(s,a) \). Show that \( J(\pi; R_\theta) = \frac{1}{1-\gamma} \boldsymbol{\theta}^\top \boldsymbol{\mu}(\pi) \), where \( \boldsymbol{\mu}(\pi) \) is the feature expectation of \( \pi \).

Solution. Plugging the linear form into the occupancy representation,

\[ \begin{aligned} J(\pi; R_\theta) &= \frac{1}{1-\gamma} \sum_{s,a} d_\pi(s,a)\, \boldsymbol{\theta}^\top \boldsymbol{\phi}(s,a) \\ &= \frac{1}{1-\gamma} \boldsymbol{\theta}^\top \sum_{s,a} d_\pi(s,a)\, \boldsymbol{\phi}(s,a) = \frac{1}{1-\gamma} \boldsymbol{\theta}^\top \boldsymbol{\mu}(\pi), \end{aligned} \]

where we swapped inner products and sums, and used the definition of \( \boldsymbol{\mu}(\pi) \).

Problem 3 (Gradient in MaxEnt IRL). Starting from the trajectory distribution \( p_{\boldsymbol{\theta}}(\tau) \propto \exp\bigl(\boldsymbol{\theta}^\top \mathbf{f}(\tau)\bigr) \), derive the gradient of the log-likelihood and show that it equals \( \hat{\boldsymbol{\mu}}_E - \boldsymbol{\mu}_{\boldsymbol{\theta}} \).

Solution. This is precisely the derivation in Section 5. Express the log-likelihood as

\[ \mathcal{L}(\boldsymbol{\theta}) = \sum_{i=1}^N \boldsymbol{\theta}^\top \mathbf{f}(\tau^{(i)}) - N \log Z(\boldsymbol{\theta}), \]

differentiate with respect to \( \boldsymbol{\theta} \), and note that \( \nabla_{\boldsymbol{\theta}} \log Z(\boldsymbol{\theta}) = \mathbb{E}_{p_{\boldsymbol{\theta}}}[\mathbf{f}(\tau)] \). Collecting terms yields the difference between empirical and model feature expectations.

Problem 4 (Max-Margin IRL as Convex Program). Consider the max-margin IRL formulation

\[ \begin{aligned} &\max_{\boldsymbol{\theta}, \delta} \;\; \delta \\ \text{s.t.}\quad &\boldsymbol{\theta}^\top \bigl( \hat{\boldsymbol{\mu}}_E - \boldsymbol{\mu}(\pi_j) \bigr) \ge \delta, \quad j=1,\dots,K, \\ &\|\boldsymbol{\theta}\|_2 \le 1. \end{aligned} \]

Prove that this is a convex optimization problem.

Solution. The feasible set is defined by linear inequalities \( \boldsymbol{\theta}^\top(\hat{\boldsymbol{\mu}}_E - \boldsymbol{\mu}(\pi_j)) - \delta \ge 0 \), which are convex, and the norm ball \( \|\boldsymbol{\theta}\|_2 \le 1 \), which is also convex. The objective \( \delta \) is linear. Therefore, we are maximizing a linear function over a convex set, i.e., a convex optimization problem.

Problem 5 (Soft Bellman Operator as Contraction). Define the soft Bellman operator \( \mathcal{T} \) on value functions by

\[ \bigl(\mathcal{T} V\bigr)(s) = \log \sum_{a \in \mathcal{A}} \exp\Bigl( R(s,a) + \gamma \sum_{s'} P(s' \mid s,a) V(s') \Bigr). \]

Show that \( \mathcal{T} \) is a contraction in the supremum norm with modulus \( \gamma \), hence soft value iteration converges to a unique fixed point.

Solution sketch. One can show that for any two value functions \( V_1, V_2 \),

\[ \bigl| (\mathcal{T} V_1)(s) - (\mathcal{T} V_2)(s) \bigr| \le \gamma \max_{s'} \bigl| V_1(s') - V_2(s') \bigr| \]

for all states \( s \), using the fact that the log-sum-exp function is 1-Lipschitz with respect to the supremum norm. Taking the supremum over \( s \) yields \( \|\mathcal{T} V_1 - \mathcal{T} V_2\|_\infty \le \gamma \|V_1 - V_2\|_\infty \). By the Banach fixed-point theorem, a unique soft value function exists and soft value iteration converges to it.

In this lesson we formalized inverse reinforcement learning for robotics. Starting from MDP preliminaries and occupancy measures, we showed how expert demonstrations define feature expectations that constrain the unknown reward function. We derived feature expectation matching and max-margin IRL, as well as maximum entropy IRL with its log-linear trajectory distribution and gradient structure. We also sketched IRL implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica, emphasizing how these integrate with standard robotics libraries and planners. In subsequent lessons and chapters we will build on these ideas to construct full imitation and reinforcement learning pipelines for advanced robotic manipulation and motion planning.