qe.solve_discrete_riccati versus scipy.linalg.solve_discrete_are

[oyamad]

Related to #84 and #356, but open a separate issue.

The claim in the docstring in solve_discrete_riccati

Note that SciPy also has a discrete riccati equation solver. However it
cannot handle the case where :math:R is not invertible, or when :math:N
is nonzero.

is no longer true.

• See this for nonzero N, and
• see below for singular R.

dare_test_tjm_2:

A = [[0, -1],
     [0, 2]]
B = [[1, 0],
     [1, 1]]
Q = [[1, 0],
     [0, 0]]
R = [[4, 2],
     [2, 1]]
A, B, Q, R = np.atleast_2d(A, B, Q, R)
# X = qe.solve_discrete_riccati(A, B, Q, R)  # LinAlgError on my osx
X_sp = scipy.linalg.solve_discrete_are(A, B, Q, R)
Y = np.zeros((2, 2))
Y[0, 0] = 1
np.max(np.abs(X_sp - Y))

1.2591976102547471e-16

dare_test_tjm_3:

r = 0.5
I = np.identity(2)
A = [[2 + r**2, 0],
     [0,        0]]
A = np.array(A)
B = I
R = [[1, r],
     [r, r*r]]
Q = I - np.dot(A.T, A) + np.dot(A.T, np.linalg.solve(R + I, A))
# X = qe.solve_discrete_riccati(A, B, Q, R)  # LinAlgError on my osx
X_sp = scipy.linalg.solve_discrete_are(A, B, Q, R)
Y = np.identity(2)
np.max(np.abs(X_sp - Y))

1.2510407154664449e-08

[ilayn]

Hi @oyamad . I was searching for something else about Riccatis but I don't know how I ended up here :) Coincidentally I'm the one who revised the SciPy version of these solvers. Please let us know if you find discrepancies in Scipy versions too. Also if you want to extend your Riccati tests, we are using the CAREX and DAREX gallery for continuous and discrete solvers which you can find here. If you follow the references you can also get the matlab mat files for testing in other domains.

By the way, very interesting project, as a control guy, I was aware of the usage of Kalman filters but I wouldn't expect to find LQ solvers inside economics package 😃 I'll definitely check it out.

[oyamad]

@ilayn Thanks for your message! We have a DARE solver based on a doubling algorithm, written by @jstac and @thomassargent30 (I am just a user, having no detailed knowledge about solution algorithms...). Our solver had a small bug, which has been fixed, and now, at least for my usecases, its outputs coincide with those by your current solve_discrete_are (up to small approximation errors).

It will be interesting to see how our solver performs for the testcases in SciPy. Thanks for the information!

[ilayn]

Ah OK. There are some generalized Lyapunov solvers planned in the pipeline. I have the direct solvers coded in some other library in case you need a basis to built on (they are 5x slow compared to matlab now as they are currently in pure Python). I will try to polish them and propose to push them into SciPy at some point. Current main problem is detecting no solution cases properly but since as far as I can see you guys are pretty into solvers maybe you might have a better idea to vectorize or increase precision.

The algorithm is the typical QZ-based symmetry exploitation and solve small linear systems of max 2x2 block size which I think what SLICOT does with extras. The solvers are in https://github.com/ilayn/harold/blob/master/harold/_solvers.py .

I don't know if MIT license is free enough but please feel free to use them at your will which was the main idea of MIT in the first place.