Feature/Sylvester solver

[stekien]

Summary

Adds sylvester_solver(A, B, C), which solves the Sylvester equation $AX + XB = C$ for $X$. Addresses #3396 (motivated by this Discourse thread).

A is n x n, B is m x m and C is n x m. Each argument can be real or complex, in any combination. The result is always a complex n x m matrix.

Primitive (stan/math/prim/fun/sylvester_solver.hpp). Bartels–Stewart with complex Schur decompositions $A = USU^H$ and $B = VTV^H$ (Eigen::ComplexSchur). This reduces the problem to $SY + YT = U^H C V$ with upper-triangular $S$ and $T$, which is solved one column at a time by triangular back substitution. Then $X = UYV^H$. The cost is $O(n^3 + m^3 + nm(n + m))$. The scalar type is templated, so forward mode (fvar) goes through this implementation.

Reverse mode (stan/math/rev/fun/sylvester_solver.hpp). The forward pass runs the same algorithm on values and keeps the Schur factors in the arena. With $G$ the adjoint of $X$ (adjoint of $\mathrm{Re},X$ plus $i$ times the adjoint of $\mathrm{Im},X$), the reverse pass solves the adjoint equation $A^H Z + Z B^H = G$ and accumulates

$$\nabla_A = -ZX^H, \qquad \nabla_B = -X^H Z, \qquad \nabla_C = Z.$$

The derivation is in the issue thread (#3396). The adjoint equation reuses the forward-pass Schur factors: it becomes $S^H W + W T^H = U^H G V$, solved column by column from last to first, so no second decomposition is needed. Complex inputs receive the real and imaginary parts of their gradient, and real inputs receive the real part. Inputs can be Matrix<var>, var_value<Matrix> or Matrix<std::complex<var>>, mixed freely with double inputs. The returned Matrix<std::complex<var>> is copied out of the arena, so writing into the result does not change what the reverse pass reads.

Lyapunov. The continuous Lyapunov equation $AX + XA^H = Q$ is the special case sylvester_solver(A, A.adjoint(), Q). This PR does not add a separate function for it. I can add a lyapunov_solver wrapper here or in a follow-up if that is preferred.

Questions for reviewers

• The result is complex even when all inputs are real. A real Schur (quasi-triangular) path would return a real matrix for real inputs, at the cost of handling 2 x 2 diagonal blocks. Is a complex result acceptable for a first version?
• A unique solution requires $\lambda_i(A) + \mu_j(B) \neq 0$ for all eigenvalues of A and B. This is not checked. If it fails, the result contains inf or NaN and no exception is thrown. Should the function check the diagonals of the shifted triangular systems and throw?
• Naming: sylvester_solver, or something closer to SciPy's solve_sylvester?

This adds the function to the Math library only. Exposing it in the Stan language needs a separate stanc3 PR.

Tests

test/unit/math/prim/fun/sylvester_test.cpp

• Real and complex systems with a known solution (6 x 5): checks the shape and residual, and that real inputs give an imaginary part at roundoff level.
• Mixed real and complex arguments.
• Diagonal A and B against the closed form $X_{ij} = C_{ij} / (a_i + b_j)$.
• Lyapunov equation $AX + XA^T = -Q$ with stable A and positive definite Q: the solution is symmetric and positive definite.
• Eigen expression arguments give the same result as plain matrices.
• Error handling: non-square A or B, and C with the wrong number of rows or columns, throw std::invalid_argument.
• Zero-size inputs return an empty result of the right shape.

test/unit/math/mix/fun/sylvester_test.cpp

• expect_ad with real inputs (A has a complex eigenvalue pair), complex inputs, and mixed real and complex inputs. This checks values and derivatives in reverse, forward and mixed mode against finite differences.
• Lyapunov case with B = A^T, so the autodiff variables of A enter twice.
• The size-mismatch error is checked for every autodiff type.
• var_value<Matrix> inputs give the same value and adjoints as Matrix<var> inputs. (expect_ad_matvar does not apply because the result has std::complex<var> scalars.)
• Writing into the returned matrix does not affect the reverse pass.

Side Effects

None for existing code. The PR adds two new headers and includes them in stan/math/prim/fun.hpp and stan/math/rev/fun.hpp.

Release notes

Added sylvester_solver(A, B, C), which solves the Sylvester equation AX + XB = C for real or complex matrices, with forward- and reverse-mode autodiff. The continuous Lyapunov equation AX + XA^H = Q can be solved as sylvester_solver(A, A^H, Q).

Checklist

• Copyright holder: Simons Foundation

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested