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, \dots\) starting from \(X_0 = B\), \(a_0 = A\):

\[a_j = a_{j-1} a_{j-1}\]
\[X_j = X_{j-1} + a_{j-1} X_{j-1} a_{j-1}'\]
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.