discrete_lyapunov_equation.cc bug in solution of Equation 4.2

[grantgib]

https://github.com/RobotLocomotion/drake/blob/4b98f4375f978a02182093b6dfabdcd3d499c236/math/discrete_lyapunov_equation.cc#L180

There is a possible bug in one of the steps for computing an intermediate term xbar. The factorization of the referenced equation 4.2 is missing an identity subtraction in the diagonal elements of MatrixXd lhs (line 199)

It should in fact be MatrixXd lhs(2 * (m - 2), 2 * (m - 2)); lhs << s_11(0, 0) * S_1.transpose() - Identity, s_11(1, 0) * S_1.transpose(), s_11(0, 1) * S_1.transpose(), s_11(1, 1) * S_1.transpose() - Identity;

[jwnimmer-tri]

\CC FYI @RussTedrake @sherm1 @hongkai-dai @FischerGundlach @weiqiao (as participants in #10885).

[grantgib]

dlyap_solution_derivation.pdf I've attached a correct derivation. The A,Q matrices used in the dlyap unit tests might be too sparse to illuminate the error. Comparing the corrected solution of Eq 4.2 with the original code returns the same result for all unit tests. Perhaps with more dense State and Cost matrices the error would become more clear.

Thanks,

[hongkai-dai]

@grantgib Thanks a lot for finding this bug and deriving the math! May I ask if you have a test example that would fail with the current code on master (without subtracting the identity)? What triggered you to find this math bug?

[grantgib]

@hongkai-dai I tried to generate some random examples that would show failure but so far both methods return the same result (6 significant digits).

As far as finding the bug, I was implementing my own solution for dlyap equations and was using this as reference.

[ggould-tri]

Tagging as mathematical_program / hongkai for now; feel free to re-tag and re-assign as appopriate.

[RussTedrake]

Amazing. @grantgib -- thanks for catching and reporting it! @hongkai-dai -- thanks for fixing it.