matrix_eqn

Filename: matrix_eqn.py

This files holds several functions that are used to solve matrix equations. Currently has functionality to solve:

  • Lyapunov Equations
  • Ricatti Equations
TODO: 1. See issue 47 on github repository, should add support for
Sylvester equations 2. Fix warnings from checking conditioning of matrices
quantecon.matrix_eqn.solve_discrete_lyapunov(A, B, max_it=50, method='doubling')[source]

Computes the solution to the discrete lyapunov equation

\[AXA' - X + B = 0\]

X is computed by using a doubling algorithm. In particular, we iterate to convergence on X_j with the following recursions for j = 1, 2,... starting from X_0 = B, a_0 = A:

\[ \begin{align}\begin{aligned}a_j = a_{j-1} a_{j-1}\\X_j = X_{j-1} + a_{j-1} X_{j-1} a_{j-1}'\end{aligned}\end{align} \]
Parameters:

A : array_like(float, ndim=2)

An n x n matrix as described above. We assume in order for convergence that the eigenvalues of A have moduli bounded by unity

B : array_like(float, ndim=2)

An n x n matrix as described above. We assume in order for convergence that the eigenvalues of A have moduli bounded by unity

max_it : scalar(int), optional(default=50)

The maximum number of iterations

method : string, optional(default=”doubling”)

Describes the solution method to use. If it is “doubling” then uses the doubling algorithm to solve, if it is “bartels-stewart” then it uses scipy’s implementation of the Bartels-Stewart approach.

Returns:

gamma1: array_like(float, ndim=2)

Represents the value V

quantecon.matrix_eqn.solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500)[source]

Solves the discrete-time algebraic Riccati equation

X = A’XA - (N + B’XA)’(B’XB + R)^{-1}(N + B’XA) + Q

via a modified structured doubling algorithm. An explanation of the algorithm can be found in the reference below.

Note that SciPy also has a discrete riccati equation solver. However it cannot handle the case where R is not invertible, or when N is nonzero. Both of these cases can be handled in the algorithm implemented below.

Parameters:

A : array_like(float, ndim=2)

k x k array.

B : array_like(float, ndim=2)

k x n array

Q : array_like(float, ndim=2)

k x k, should be symmetric and non-negative definite

R : array_like(float, ndim=2)

n x n, should be symmetric and positive definite

N : array_like(float, ndim=2)

n x k array

tolerance : scalar(float), optional(default=1e-10)

The tolerance level for convergence

max_iter : scalar(int), optional(default=500)

The maximum number of iterations allowed

Returns:

X : array_like(float, ndim=2)

The fixed point of the Riccati equation; a k x k array representing the approximate solution

References

Chiang, Chun-Yueh, Hung-Yuan Fan, and Wen-Wei Lin. “STRUCTURED DOUBLING ALGORITHM FOR DISCRETE-TIME ALGEBRAIC RICCATI EQUATIONS WITH SINGULAR CONTROL WEIGHTING MATRICES.” Taiwanese Journal of Mathematics 14, no. 3A (2010): pp-935.