# 63. Lagrangian for LQ Control¶

```
!pip install quantecon
```

```
Requirement already satisfied: quantecon in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (0.6.0)
Requirement already satisfied: scipy>=1.5.0 in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from quantecon) (1.9.1)
Requirement already satisfied: numba in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from quantecon) (0.56.4)
Requirement already satisfied: requests in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from quantecon) (2.28.1)
Requirement already satisfied: sympy in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from quantecon) (1.10.1)
Requirement already satisfied: numpy>=1.17.0 in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from quantecon) (1.23.5)
Requirement already satisfied: llvmlite<0.40,>=0.39.0dev0 in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from numba->quantecon) (0.39.1)
Requirement already satisfied: setuptools in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from numba->quantecon) (63.4.1)
Requirement already satisfied: idna<4,>=2.5 in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from requests->quantecon) (3.3)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from requests->quantecon) (1.26.11)
Requirement already satisfied: certifi>=2017.4.17 in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from requests->quantecon) (2022.9.14)
Requirement already satisfied: charset-normalizer<3,>=2 in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from requests->quantecon) (2.0.4)
Requirement already satisfied: mpmath>=0.19 in /__w/lecture-python.myst/lecture-python.myst/3/envs/quantecon/lib/python3.9/site-packages (from sympy->quantecon) (1.2.1)
```

```
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv
```

```
import numpy as np
from quantecon import LQ
from scipy.linalg import schur
```

## 63.1. Overview¶

This is a sequel to this lecture linear quadratic dynamic programming

It can also be regarded as presenting **invariant subspace** techniques that extend ones that
we encountered earlier in this lecture stability in linear rational expectations models

We present a Lagrangian formulation of an infinite horizon linear quadratic undiscounted dynamic programming problem.

Such a problem is also sometimes called an optimal linear regulator problem.

A Lagrangian formulation

carries insights about connections between stability and optimality

is the basis for fast algorithms for solving Riccati equations

opens the way to constructing solutions of dynamic systems that don’t come directly from an intertemporal optimization problem

A key tool in this lecture is the concept of an \(n \times n\) **symplectic** matrix.

A symplectic matrix has eigenvalues that occur in **reciprocal pairs**, meaning that if \(\lambda_i \in (-1,1)\) is an eigenvalue, then so is \(\lambda_i^{-1}\).

This reciprocal pairs property of the eigenvalues of a matrix is a tell-tale sign that the matrix describes the joint dynamics of a system of
equations describing the **states** and **costates** that constitute first-order necessary conditions for solving an undiscounted linear-quadratic infinite-horizon optimization problem.

The symplectic matrix that will interest us describes the first-order dynamics of **state** and **co-state** vectors of an optimally controlled system.

In focusing on eigenvalues and eigenvectors of this matrix, we capitalize on an analysis of
**invariant subspaces.**

These invariant subspace formulations of LQ dynamic programming problems provide a bridge between recursive (i.e., dynamic programming) formulations and classical formulations of linear control and linear filtering problems that make use of related matrix decompositions (see for example this lecture and this lecture).

While most of this lecture focuses on undiscounted problems, later sections describe handy ways of transforming discounted problems to undiscounted ones.

The techniques in this lecture will prove useful when we study Stackelberg and Ramsey problem in this lecture.

## 63.2. Undiscounted LQ DP Problem¶

The problem is to choose a sequence of controls \(\{u_t\}_{t=0}^\infty\) to maximize the criterion

subject to \(x_{t+1}=Ax_t+Bu_t\), where \(x_0\) is a given initial state vector.

Here \(x_t\) is an \((n\times 1)\) vector of state variables, \(u_t\) is a \((k\times 1)\) vector of controls, \(R\) is a positive semidefinite symmetric matrix, \(Q\) is a positive definite symmetric matrix, \(A\) is an \((n\times n)\) matrix, and \(B\) is an \((n\times k)\) matrix.

The optimal value function turns out to be quadratic, \(V(x)= - x'Px\), where \(P\) is a positive semidefinite symmetric matrix.

Using the transition law to eliminate next period’s state, the Bellman equation becomes

The first-order necessary conditions for the maximum problem on the right side of equation (63.1) are

Note

We use the following rules for differentiating quadratic and bilinear matrix forms: \({\partial x' A x \over \partial x} = (A + A') x; {\partial y' B z \over \partial y} = B z, {\partial y' B z \over \partial z} = B' y\).

which implies that an optimal decision rule for \(u\) is

or

where

Substituting \(u = - (Q+B'PB)^{-1}B'PAx\) into the right side of equation (63.1) and rearranging gives

Equation (63.2) is called an **algebraic matrix Riccati** equation.

There are multiple solutions of equation (63.2).

But only one of them is positive definite.

The positive define solution is associated with the maximum of our problem.

It expresses the matrix \(P\) as an implicit function of the matrices \(R,Q,A,B\).

Notice that the **gradient of the value function** is

We shall use fact (63.3) later.

## 63.3. Lagrangian¶

For the undiscounted optimal linear regulator problem, form the Lagrangian

where \(2 \mu_{t+1}\) is a vector of Lagrange multipliers on the time \(t\) transition law \(x_{t+1} = A x_t + B u_t\).

(We put the \(2\) in front of \(\mu_{t+1}\) to make things match up nicely with equation (63.3).)

First-order conditions for maximization with respect to \(\{u_t,x_{t+1}\}_{t=0}^\infty\) are

Define \(\mu_0\) to be a vector of shadow prices of \(x_0\) and apply an envelope condition to (63.4) to deduce that

which is a time \(t=0 \) counterpart to the second equation of system (63.5).

An important fact is that

where \(P\) is a positive define matrix that solves the algebraic Riccati equation (63.2).

Thus, from equations (63.3) and (63.6), \(- 2 \mu_{t}\) is the gradient of the value function with respect to \(x_t\).

The Lagrange multiplier vector \(\mu_{t}\) is often called the **costate** vector that
corresponds to the **state** vector \(x_t\).

It is useful to proceed with the following steps:

solve the first equation of (63.5) for \(u_t\) in terms of \(\mu_{t+1}\).

substitute the result into the law of motion \(x_{t+1} = A x_t + B u_t\).

arrange the resulting equation and the second equation of (63.5) into the form

where

When \(L\) is of full rank (i.e., when \(A\) is of full rank), we can write system (63.7) as

where

## 63.4. State-Costate Dynamics¶

We seek to solve the difference equation system (63.8) for a sequence \(\{x_t\}_{t=0}^\infty\) that satisfies

an initial condition for \(x_0\)

a terminal condition \(\lim_{t \rightarrow +\infty} x_t =0\)

This terminal condition reflects our desire for a **stable** solution, one that does not diverge as \(t \rightarrow \infty\).

We inherit our wish for stability of the \(\{x_t\}\) sequence from a desire to maximize

which requires that \(x_t' R x_t\) converge to zero as \(t \rightarrow + \infty\).

## 63.5. Reciprocal Pairs Property¶

To proceed, we study properties of the \((2n \times 2n)\) matrix \(M\) defined in (63.9).

It helps to introduce a \((2n \times 2n)\) matrix

The rank of \(J\) is \(2n\).

**Definition:** A matrix \(M\) is called **symplectic** if

Salient properties of symplectic matrices that are readily verified include:

If \(M\) is symplectic, then \(M^2\) is symplectic

The determinant of a symplectic, then \(\textrm{det}(M) = 1\)

It can be verified directly that \(M\) in equation (63.9) is symplectic.

It follows from equation (63.10) and from the fact \(J^{-1} = J^\prime = -J\) that for any symplectic matrix \(M\),

Equation (63.11) states that \(M^\prime\) is related to the inverse of \(M\)
by a **similarity transformation**.

For square matrices, recall that

similar matrices share eigenvalues

eigenvalues of the inverse of a matrix are inverses of eigenvalues of the matrix

a matrix and its transpose share eigenvalues

It then follows from equation (63.11) that the eigenvalues of \(M\) occur in reciprocal pairs: if \(\lambda\) is an eigenvalue of \(M\), so is \(\lambda^{-1}\).

Write equation (63.8) as

where \(y_t = \begin{pmatrix}x_t\cr \mu_t\cr\end{pmatrix}\).

Consider a **triangularization** of \(M\)

where

each block on the right side is \((n\times n)\)

\(V\) is nonsingular

all eigenvalues of \(W_{22}\) exceed \(1\) in modulus

all eigenvalues of \(W_{11}\) are less than \(1\) in modulus

## 63.6. Schur decomposition¶

The **Schur decomposition** and the **eigenvalue decomposition**
are two decompositions of the form (63.13).

Write equation (63.12) as

A solution of equation (63.14) for arbitrary initial condition \(y_0\) is evidently

where \(W_{12,t} = W_{12}\) for \(t=1\) and for \(t \geq 2\) obeys the recursion

and where \(W^t_{ii}\) is \(W_{ii}\) raised to the \(t\)th power.

Write equation (63.15) as

where \(y^\ast_t = V^{-1} y_t\), and in particular where

and where \(V^{ij}\) denotes the \((i,j)\) piece of the partitioned \(V^{-1}\) matrix.

Because \(W_{22}\) is an unstable matrix, \(y^\ast_t\) will diverge unless \(y^\ast_{20} = 0\).

Let \(V^{ij}\) denote the \((i,j)\) piece of the partitioned \(V^{-1}\) matrix.

To attain stability, we must impose \(y^\ast_{20} =0\), which from equation (63.16) implies

or

This equation replicates itself over time in the sense that it implies

But notice that because \((V^{21}\ V^{22})\) is the second row block of the inverse of \(V,\) it follows that

which implies

Therefore,

So we can write

and

However, we know that \(\mu_t = P x_t\), where \(P\) occurs in the matrix that solves the Riccati equation.

Thus, the preceding argument establishes that

Remarkably, formula (63.17) provides us with a computationally efficient way of computing the positive definite matrix \(P\) that solves the algebraic Riccati equation (63.2) that emerges from dynamic programming.

This same method can be applied to compute the solution of any system of the form (63.8) if a solution exists, even if eigenvalues of \(M\) fail to occur in reciprocal pairs.

The method will typically work so long as the eigenvalues of \(M\) split half inside and half outside the unit circle.

Systems in which eigenvalues (properly adjusted for discounting) fail to occur in reciprocal pairs arise when the system being solved is an equilibrium of a model in which there are distortions that prevent there being any optimum problem that the equilibrium solves. See [LS18], ch 12.

## 63.7. Application¶

Here we demonstrate the computation with an example which is the deterministic version of an example borrowed from this quantecon lecture.

```
# Model parameters
r = 0.05
c_bar = 2
μ = 1
# Formulate as an LQ problem
Q = np.array([[1]])
R = np.zeros((2, 2))
A = [[1 + r, -c_bar + μ],
[0, 1]]
B = [[-1],
[0]]
# Construct an LQ instance
lq = LQ(Q, R, A, B)
```

Given matrices \(A\), \(B\), \(Q\), \(R\), we can then compute \(L\), \(N\), and \(M=L^{-1}N\).

```
def construct_LNM(A, B, Q, R):
n, k = lq.n, lq.k
# construct L and N
L = np.zeros((2*n, 2*n))
L[:n, :n] = np.eye(n)
L[:n, n:] = B @ np.linalg.inv(Q) @ B.T
L[n:, n:] = A.T
N = np.zeros((2*n, 2*n))
N[:n, :n] = A
N[n:, :n] = -R
N[n:, n:] = np.eye(n)
# compute M
M = np.linalg.inv(L) @ N
return L, N, M
```

```
L, N, M = construct_LNM(lq.A, lq.B, lq.Q, lq.R)
```

```
M
```

```
array([[ 1.05 , -1. , -0.95238095, 0. ],
[ 0. , 1. , 0. , 0. ],
[ 0. , 0. , 0.95238095, 0. ],
[ 0. , 0. , 0.95238095, 1. ]])
```

Let’s verify that \(M\) is symplectic.

```
n = lq.n
J = np.zeros((2*n, 2*n))
J[n:, :n] = np.eye(n)
J[:n, n:] = -np.eye(n)
M @ J @ M.T - J
```

```
array([[-1.32169408e-17, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00]])
```

We can compute the eigenvalues of \(M\) using `np.linalg.eigvals`

, arranged in ascending order.

```
eigvals = sorted(np.linalg.eigvals(M))
eigvals
```

```
[0.9523809523809523, 1.0, 1.0, 1.05]
```

When we apply Schur decomposition such that \(M=V W V^{-1}\), we want

the upper left block of \(W\), \(W_{11}\), to have all of its eigenvalues less than 1 in modulus, and

the lower right block \(W_{22}\) to have eigenvalues that exceed 1 in modulus.

To get what we want, let’s define a sorting function that tells `scipy.schur`

to sort the corresponding eigenvalues with modulus smaller than 1 to the upper left.

```
stable_eigvals = eigvals[:n]
def sort_fun(x):
"Sort the eigenvalues with modules smaller than 1 to the top-left."
if x in stable_eigvals:
stable_eigvals.pop(stable_eigvals.index(x))
return True
else:
return False
W, V, _ = schur(M, sort=sort_fun)
```

```
W
```

```
array([[ 1. , -0.02316402, -1.00085948, -0.95000594],
[ 0. , 0.95238095, -0.00237501, -0.95325452],
[ 0. , 0. , 1.05 , 0.02432222],
[ 0. , 0. , 0. , 1. ]])
```

```
V
```

```
array([[ 0.99875234, 0.00121459, -0.04992284, 0. ],
[ 0.04993762, -0.02429188, 0.99845688, 0. ],
[ 0. , 0.04992284, 0.00121459, 0.99875234],
[ 0. , -0.99845688, -0.02429188, 0.04993762]])
```

We can check the modulus of eigenvalues of \(W_{11}\) and \(W_{22}\).

Since they are both triangular matrices, eigenvalues are the diagonal elements.

```
# W11
np.diag(W[:n, :n])
```

```
array([1. , 0.95238095])
```

```
# W22
np.diag(W[n:, n:])
```

```
array([1.05, 1. ])
```

The following functions wrap \(M\) matrix construction, Schur decomposition, and stability-imposing computation of \(P\).

```
def stable_solution(M, verbose=True):
"""
Given a system of linear difference equations
y' = |a b| y
x' = |c d| x
which is potentially unstable, find the solution
by imposing stability.
Parameter
---------
M : np.ndarray(float)
The matrix represents the linear difference equations system.
"""
n = M.shape[0] // 2
stable_eigvals = list(sorted(np.linalg.eigvals(M))[:n])
def sort_fun(x):
"Sort the eigenvalues with modules smaller than 1 to the top-left."
if x in stable_eigvals:
stable_eigvals.pop(stable_eigvals.index(x))
return True
else:
return False
W, V, _ = schur(M, sort=sort_fun)
if verbose:
print('eigenvalues:\n')
print(' W11: {}'.format(np.diag(W[:n, :n])))
print(' W22: {}'.format(np.diag(W[n:, n:])))
# compute V21 V11^{-1}
P = V[n:, :n] @ np.linalg.inv(V[:n, :n])
return W, V, P
def stationary_P(lq, verbose=True):
"""
Computes the matrix :math:`P` that represent the value function
V(x) = x' P x
in the infinite horizon case. Computation is via imposing stability
on the solution path and using Schur decomposition.
Parameters
----------
lq : qe.LQ
QuantEcon class for analyzing linear quadratic optimal control
problems of infinite horizon form.
Returns
-------
P : array_like(float)
P matrix in the value function representation.
"""
Q = lq.Q
R = lq.R
A = lq.A * lq.beta ** (1/2)
B = lq.B * lq.beta ** (1/2)
n, k = lq.n, lq.k
L, N, M = construct_LNM(A, B, Q, R)
W, V, P = stable_solution(M, verbose=verbose)
return P
```

```
# compute P
stationary_P(lq)
```

```
eigenvalues:
W11: [1. 0.95238095]
W22: [1.05 1. ]
```

```
array([[ 0.1025, -2.05 ],
[-2.05 , 41. ]])
```

Note that the matrix \(P\) computed in this way is close to what we get from the routine in quantecon that solves an algebraic Riccati equation by iterating to convergence on a Riccati difference equation.

The small difference comes from computational errors and will decrease as we increase the maximum number of iterations or decrease the tolerance for convergence.

```
lq.stationary_values()
```

```
(array([[ 0.1025, -2.05 ],
[-2.05 , 41.01 ]]),
array([[-0.09761905, 1.95238095]]),
0)
```

Using a Schur decomposition is much more efficient.

```
%%timeit
stationary_P(lq, verbose=False)
```

```
169 µs ± 347 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
```

```
%%timeit
lq.stationary_values()
```

```
2.89 ms ± 8.84 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```

## 63.8. Other Applications¶

The preceding approach to imposing stability on a system of potentially unstable linear difference equations is not limited to linear quadratic dynamic optimization problems.

For example, the same method is used in our Stability in Linear Rational Expectations Models lecture.

Let’s try to solve the model described in that lecture by applying the `stable_solution`

function defined in this lecture above.

```
def construct_H(ρ, λ, δ):
"contruct matrix H given parameters."
H = np.empty((2, 2))
H[0, :] = ρ,δ
H[1, :] = - (1 - λ) / λ, 1 / λ
return H
H = construct_H(ρ=.9, λ=.5, δ=0)
```

```
W, V, P = stable_solution(H)
P
```

```
eigenvalues:
W11: [0.9]
W22: [2.]
```

```
array([[0.90909091]])
```

## 63.9. Discounted Problems¶

### 63.9.1. Transforming States and Controls to Eliminate Discounting¶

A pair of useful transformations allows us to convert a discounted problem into an undiscounted one.

Thus, suppose that we have a discounted problem with objective

and that the state transition equation is again \(x_{t +1 }=Ax_t+Bu_t\).

Define the transformed state and control variables

\(\hat x_t = \beta^{\frac{t}{2}} x_t \)

\(\hat u_t = \beta^{\frac{t}{2}} u_t\)

and the transformed transition equation matrices

\(\hat A = \beta^{\frac{1}{2}} A\)

\(\hat B = \beta^{\frac{1}{2}} B \)

so that the adjusted state and control variables obey the transition law

Then a discounted optimal control problem defined by \(A, B, R, Q, \beta\) having optimal policy characterized by \(P, F\) is associated with an equivalent undiscounted problem defined by \(\hat A, \hat B, Q, R\) having optimal policy characterized by \(\hat F, \hat P\) that satisfy the following equations:

and

It follows immediately from the definitions of \(\hat A, \hat B\) that \(\hat F = F\) and \(\hat P = P\).

By exploiting these transformations, we can solve a discounted problem by solving an associated undiscounted problem.

In particular, we can first transform a discounted LQ problem to an undiscounted one and then solve that discounted optimal regulator problem using the Lagrangian and invariant subspace methods described above.

For example, when \(\beta=\frac{1}{1+r}\), we can solve for \(P\) with \(\hat{A}=\beta^{1/2} A\) and \(\hat{B}=\beta^{1/2} B\).

These settings are adopted by default in the function `stationary_P`

defined above.

```
β = 1 / (1 + r)
lq.beta = β
```

```
stationary_P(lq)
```

```
eigenvalues:
W11: [0.97590007 0.97590007]
W22: [1.02469508 1.02469508]
```

```
array([[ 0.0525, -1.05 ],
[-1.05 , 21. ]])
```

We can verify that the solution agrees with one that comes from applying the routine `LQ.stationary_values`

in the quantecon package.

```
lq.stationary_values()
```

```
(array([[ 0.0525, -1.05 ],
[-1.05 , 21. ]]),
array([[-0.05, 1. ]]),
0.0)
```

### 63.9.2. Lagrangian for Discounted Problem¶

For several purposes, it is useful explicitly briefly to describe a Lagrangian for a discounted problem.

Thus, for the discounted optimal linear regulator problem, form the Lagrangian

where \(2 \mu_{t+1}\) is a vector of Lagrange multipliers on the state vector \(x_{t+1}\).

First-order conditions for maximization with respect to \(\{u_t,x_{t+1}\}_{t=0}^\infty\) are

Define \(2 \mu_0\) to be the vector of shadow prices of \(x_0\) and apply an envelope condition to (63.18) to deduce that

which is a time \(t=0 \) counterpart to the second equation of system (63.19).

Proceeding as we did above with the undiscounted system (63.5), we can rearrange the first-order conditions into the system

which in the special case that \(\beta = 1\) agrees with equation (63.5), as expected.

By staring at system (63.20), we can infer identities that shed light on the structure of optimal linear regulator problems, some of which will be useful in this lecture when we apply and extend the methods of this lecture to study Stackelberg and Ramsey problems.

First, note that the first block of equation system (63.20) asserts that when \(\mu_{t+1} = P x_{t+1}\), then

which can be rearranged to sbe

This expression for the optimal closed loop dynamics of the state must agree with an alternative expression that we had derived with dynamic programming, namely,

But using

it follows that

Thus, our two expressions for the closed loop dynamics agree if and only if

Matrix equation (63.22) can be verified by applying a partitioned inverse formula.

Note

Just use the formula \((a - b d^{-1} c)^{-1} = a^{-1} + a^{-1} b (d - c a^{-1} b)^{-1} c a^{-1}\) for appropriate choices of the matrices \(a, b, c, d\).

Next, note that for *any* fixed \(F\) for which eigenvalues of \(A- BF\) are less than \(\frac{1}{\beta}\) in modulus, the value function associated with using this rule forever is \(- x_0 \tilde P x_0\) where \(\tilde P\) obeys the following matrix equation:

Evidently, \(\tilde P = P \) only when \(F \) obeys formula (63.21).

Next, note that the second equation of system (63.20) implies the “forward looking” equation for the Lagrange multiplier

whose solution is

where

where we must require that \(F\) obeys equation (63.21).

Equations (63.23) and (63.24) provide different perspectives on the optimal value function.