From c136695e2800a4b316bc359953aa2ea7840c846d Mon Sep 17 00:00:00 2001 From: Stefan Kienle <127690239+stekien@users.noreply.github.com> Date: Thu, 15 Oct 2026 14:34:36 -0400 Subject: [PATCH 1/6] add sylvester solver (AX + XB = C) and tests --- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/sylvester_solver.hpp | 72 ++++++++++ test/unit/math/prim/fun/sylvester_test.cpp | 158 +++++++++++++++++++++ 3 files changed, 231 insertions(+) create mode 100644 stan/math/prim/fun/sylvester_solver.hpp create mode 100644 test/unit/math/prim/fun/sylvester_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 8ef4204c81b..cdd963d31d4 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -313,6 +313,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/sylvester_solver.hpp b/stan/math/prim/fun/sylvester_solver.hpp new file mode 100644 index 00000000000..179a6048bfe --- /dev/null +++ b/stan/math/prim/fun/sylvester_solver.hpp @@ -0,0 +1,72 @@ +#ifndef STAN_MATH_PRIM_FUN_SYLVESTER_SOLVER_HPP +#define STAN_MATH_PRIM_FUN_SYLVESTER_SOLVER_HPP + +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Solve the Sylvester equation A X + X B = C for X. + * + * With complex Schur decompositions A = U S U^H and B = V T V^H the + * equation reduces to S Y + Y T = U^H C V, solved column by column, + * and X = U Y V^H. A unique solution requires lambda_i(A) + mu_j(B) != 0 + * for all eigenvalues of A and B. + * + * @tparam MatA type of A (real or complex scalars) + * @tparam MatB type of B (real or complex scalars) + * @tparam MatC type of C (real or complex scalars) + * @param[in] A n x n matrix + * @param[in] B m x m matrix + * @param[in] C n x m matrix + * @return complex n x m solution X + * @throw std::invalid_argument if A or B is not square, or if the + * dimensions of C do not match + */ +template * = nullptr, + require_all_t>, + std::is_arithmetic>, + std::is_arithmetic>>* = nullptr> +inline auto sylvester_solver(const MatA& A, const MatB& B, + const MatC& C) { + check_square("sylvester_solver", "A", A); + check_square("sylvester_solver", "B", B); + check_size_match("sylvester_solver", "rows of A", A.rows(), "rows of C", + C.rows()); + check_size_match("sylvester_solver", "rows of B", B.rows(), "cols of C", + C.cols()); + + const Eigen::Index m = B.rows(); + + Eigen::ComplexSchur> schur_A(A); + const auto& U = schur_A.matrixU(); + const auto& S = schur_A.matrixT(); + + Eigen::ComplexSchur> schur_B(B); + const auto& V = schur_B.matrixU(); + const auto& T = schur_B.matrixT(); + + // Y holds D = U^H C V initially, overwritten column by column + auto Y = (U.adjoint() * C * V).eval(); + auto Sshift = S; + for (Eigen::Index i = 0; i < m; ++i) { + Y.col(i).noalias() -= Y.leftCols(i) * T.col(i).head(i); + Sshift.diagonal() = S.diagonal().array() + T(i, i); + Sshift.template triangularView().solveInPlace(Y.col(i)); + } + + auto X = (U * Y * V.adjoint()).eval(); + + return X; +} + +} // namespace math +} // namespace stan +#endif \ No newline at end of file diff --git a/test/unit/math/prim/fun/sylvester_test.cpp b/test/unit/math/prim/fun/sylvester_test.cpp new file mode 100644 index 00000000000..2e49bd24386 --- /dev/null +++ b/test/unit/math/prim/fun/sylvester_test.cpp @@ -0,0 +1,158 @@ +#include +#include +#include +#include + +// Diagonal shifts keep lambda_i(A) + mu_j(B) well away from zero, so the +// problems are well conditioned. Random nonsymmetric matrices still have +// complex eigenvalue pairs, so the complex Schur path is exercised. + +TEST(MathFunctions, sylvester_solver_real_known_solution) { + using stan::math::sylvester_solver; + const int n = 6; + const int m = 5; + Eigen::MatrixXd A = Eigen::MatrixXd::Random(n, n) + + 2.0 * n * Eigen::MatrixXd::Identity(n, n); + Eigen::MatrixXd B = Eigen::MatrixXd::Random(m, m) + + 2.0 * m * Eigen::MatrixXd::Identity(m, m); + Eigen::MatrixXd X_true = Eigen::MatrixXd::Random(n, m); + Eigen::MatrixXd C = A * X_true + X_true * B; + + Eigen::MatrixXcd X = sylvester_solver(A, B, C); + + EXPECT_EQ(X.rows(), n); + EXPECT_EQ(X.cols(), m); + // real inputs give a real solution up to roundoff + EXPECT_LT(X.imag().norm(), 1e-12); + EXPECT_LT((X.real() - X_true).norm(), 1e-10); + EXPECT_LT((A * X + X * B - C).norm(), 1e-10); +} + +TEST(MathFunctions, sylvester_solver_complex_known_solution) { + using stan::math::sylvester_solver; + const int n = 6; + const int m = 5; + Eigen::MatrixXcd A = Eigen::MatrixXcd::Random(n, n) + + 2.0 * n * Eigen::MatrixXcd::Identity(n, n); + Eigen::MatrixXcd B = Eigen::MatrixXcd::Random(m, m) + + 2.0 * m * Eigen::MatrixXcd::Identity(m, m); + Eigen::MatrixXcd X_true = Eigen::MatrixXcd::Random(n, m); + Eigen::MatrixXcd C = A * X_true + X_true * B; + + Eigen::MatrixXcd X = sylvester_solver(A, B, C); + + EXPECT_LT((X - X_true).norm(), 1e-10); + EXPECT_LT((A * X + X * B - C).norm(), 1e-10); +} + +TEST(MathFunctions, sylvester_solver_mixed_real_complex) { + using stan::math::sylvester_solver; + const int n = 4; + const int m = 3; + Eigen::MatrixXd A_r = Eigen::MatrixXd::Random(n, n) + + 2.0 * n * Eigen::MatrixXd::Identity(n, n); + Eigen::MatrixXd B_r = Eigen::MatrixXd::Random(m, m) + + 2.0 * m * Eigen::MatrixXd::Identity(m, m); + Eigen::MatrixXcd A_c = Eigen::MatrixXcd::Random(n, n) + + 2.0 * n * Eigen::MatrixXcd::Identity(n, n); + + // real A and B, complex C: known solution + Eigen::MatrixXcd X_true = Eigen::MatrixXcd::Random(n, m); + Eigen::MatrixXcd C_c = A_r * X_true + X_true * B_r; + Eigen::MatrixXcd X1 = sylvester_solver(A_r, B_r, C_c); + EXPECT_LT((X1 - X_true).norm(), 1e-10); + + // complex A, real B and C: residual + Eigen::MatrixXd C_r = Eigen::MatrixXd::Random(n, m); + Eigen::MatrixXcd X2 = sylvester_solver(A_c, B_r, C_r); + EXPECT_LT((A_c * X2 + X2 * B_r - C_r).norm(), 1e-10); +} + +TEST(MathFunctions, sylvester_solver_diagonal) { + // for diagonal A and B: X(i, j) = C(i, j) / (A(i, i) + B(j, j)) + Eigen::VectorXd a(3); + a << 1, 2, 3; + Eigen::VectorXd b(2); + b << 4, 5; + Eigen::MatrixXd A = a.asDiagonal(); + Eigen::MatrixXd B = b.asDiagonal(); + Eigen::MatrixXd C(3, 2); + C << 1, 2, 3, 4, 5, 6; + + Eigen::MatrixXcd X = stan::math::sylvester_solver(A, B, C); + + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 2; ++j) { + EXPECT_NEAR(X(i, j).real(), C(i, j) / (a(i) + b(j)), 1e-12); + EXPECT_NEAR(X(i, j).imag(), 0.0, 1e-12); + } + } +} + +TEST(MathFunctions, sylvester_solver_lyapunov) { + // continuous Lyapunov equation A X + X A^T = -Q + using stan::math::sylvester_solver; + const int n = 5; + // stable A: by Gershgorin all eigenvalues have Re <= -(n + 1) < 0, + // so lambda_i + lambda_j != 0 and the problem is well posed + Eigen::MatrixXd A = Eigen::MatrixXd::Random(n, n) + - 2.0 * n * Eigen::MatrixXd::Identity(n, n); + Eigen::MatrixXd L = Eigen::MatrixXd::Random(n, n); + Eigen::MatrixXd Q = L * L.transpose() + Eigen::MatrixXd::Identity(n, n); + + Eigen::MatrixXcd X = sylvester_solver(A, A.transpose(), -Q); + Eigen::MatrixXd Xr = X.real(); + + EXPECT_LT(X.imag().norm(), 1e-12); + EXPECT_LT((A * Xr + Xr * A.transpose() + Q).norm(), 1e-10); + + // symmetric C gives a symmetric solution + EXPECT_LT((Xr - Xr.transpose()).norm(), 1e-12); + + // stable A and positive definite Q give a positive definite solution + Eigen::SelfAdjointEigenSolver eig( + 0.5 * (Xr + Xr.transpose())); + EXPECT_GT(eig.eigenvalues().minCoeff(), 0.0); +} + +TEST(MathFunctions, sylvester_solver_expression_inputs) { + // Eigen expressions go through plain_type_t; result must match plain inputs + using stan::math::sylvester_solver; + Eigen::MatrixXd A = Eigen::MatrixXd::Random(4, 4) + + 8.0 * Eigen::MatrixXd::Identity(4, 4); + Eigen::MatrixXd B = Eigen::MatrixXd::Random(3, 3) + + 6.0 * Eigen::MatrixXd::Identity(3, 3); + Eigen::MatrixXd C = Eigen::MatrixXd::Random(4, 3); + Eigen::MatrixXd At = A.transpose(); + + Eigen::MatrixXcd X_expr = sylvester_solver(A.transpose(), B, 2.0 * C); + Eigen::MatrixXcd X_plain = sylvester_solver(At, B, Eigen::MatrixXd(2.0 * C)); + + EXPECT_LT((X_expr - X_plain).norm(), 1e-14); +} + +TEST(MathFunctions, sylvester_solver_errors) { + using stan::math::sylvester_solver; + Eigen::MatrixXd A = Eigen::MatrixXd::Identity(3, 3); + Eigen::MatrixXd B = Eigen::MatrixXd::Identity(2, 2); + Eigen::MatrixXd C = Eigen::MatrixXd::Ones(3, 2); + + // valid input does not throw + EXPECT_NO_THROW(sylvester_solver(A, B, C)); + + // A not square + Eigen::MatrixXd A_rect = Eigen::MatrixXd::Ones(3, 2); + EXPECT_THROW(sylvester_solver(A_rect, B, C), std::invalid_argument); + + // B not square + Eigen::MatrixXd B_rect = Eigen::MatrixXd::Ones(2, 3); + EXPECT_THROW(sylvester_solver(A, B_rect, C), std::invalid_argument); + + // rows of C do not match A + Eigen::MatrixXd C_rows = Eigen::MatrixXd::Ones(4, 2); + EXPECT_THROW(sylvester_solver(A, B, C_rows), std::invalid_argument); + + // cols of C do not match B + Eigen::MatrixXd C_cols = Eigen::MatrixXd::Ones(3, 3); + EXPECT_THROW(sylvester_solver(A, B, C_cols), std::invalid_argument); +} \ No newline at end of file From 9a302d927e113d940c6443c7f598edc02ecd7102 Mon Sep 17 00:00:00 2001 From: Stefan Kienle <127690239+stekien@users.noreply.github.com> Date: Thu, 15 Oct 2026 17:10:46 -0400 Subject: [PATCH 2/6] generalize prim function to allow fvar types --- stan/math/prim/fun/sylvester_solver.hpp | 36 ++++++++++++++++++++----- 1 file changed, 30 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/fun/sylvester_solver.hpp b/stan/math/prim/fun/sylvester_solver.hpp index 179a6048bfe..ebac314509a 100644 --- a/stan/math/prim/fun/sylvester_solver.hpp +++ b/stan/math/prim/fun/sylvester_solver.hpp @@ -31,11 +31,28 @@ namespace math { */ template * = nullptr, - require_all_t>, - std::is_arithmetic>, - std::is_arithmetic>>* = nullptr> + require_all_not_t>, + is_var>, + is_var>>* = nullptr> inline auto sylvester_solver(const MatA& A, const MatB& B, const MatC& C) { + + // common base scalar: double, or fvar<...> if any input is fvar + using base_t = return_type_t, base_type_t, + base_type_t>; + using complex_t = std::complex; + using complex_mat_t + = Eigen::Matrix; + // A and B are promoted to base_t but stay real or complex as given + using A_mat_t = Eigen::Matrix< + std::conditional_t>::value, complex_t, + base_t>, + Eigen::Dynamic, Eigen::Dynamic>; + using B_mat_t = Eigen::Matrix< + std::conditional_t>::value, complex_t, + base_t>, + Eigen::Dynamic, Eigen::Dynamic>; + check_square("sylvester_solver", "A", A); check_square("sylvester_solver", "B", B); check_size_match("sylvester_solver", "rows of A", A.rows(), "rows of C", @@ -43,13 +60,20 @@ inline auto sylvester_solver(const MatA& A, const MatB& B, check_size_match("sylvester_solver", "rows of B", B.rows(), "cols of C", C.cols()); + if(unlikely(A.rows() == 0 || B.rows() == 0)) { + return complex_mat_t(A.rows(), B.rows()); + } + const Eigen::Index m = B.rows(); - Eigen::ComplexSchur> schur_A(A); + // ComplexSchur needs a plain matrix (same pattern as complex_schur_decompose) + Eigen::ComplexSchur schur_A{ + A_mat_t(A.template cast())}; const auto& U = schur_A.matrixU(); const auto& S = schur_A.matrixT(); - Eigen::ComplexSchur> schur_B(B); + Eigen::ComplexSchur schur_B{ + B_mat_t(B.template cast())}; const auto& V = schur_B.matrixU(); const auto& T = schur_B.matrixT(); @@ -62,7 +86,7 @@ inline auto sylvester_solver(const MatA& A, const MatB& B, Sshift.template triangularView().solveInPlace(Y.col(i)); } - auto X = (U * Y * V.adjoint()).eval(); + complex_mat_t X = U * Y * V.adjoint(); return X; } From 7bddb58bbf7f8d2a071540fd6e07bd9af78e61a4 Mon Sep 17 00:00:00 2001 From: Stefan Kienle <127690239+stekien@users.noreply.github.com> Date: Thu, 15 Oct 2026 17:49:20 -0400 Subject: [PATCH 3/6] add reev version of the sylvester solver and tests --- stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/sylvester_solver.hpp | 146 +++++++++++++++++++++ test/unit/math/mix/fun/sylvester_test.cpp | 147 ++++++++++++++++++++++ 3 files changed, 294 insertions(+) create mode 100644 stan/math/rev/fun/sylvester_solver.hpp create mode 100644 test/unit/math/mix/fun/sylvester_test.cpp diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index e5ef938a302..e769943e6b3 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -163,6 +163,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/sylvester_solver.hpp b/stan/math/rev/fun/sylvester_solver.hpp new file mode 100644 index 00000000000..539e7aa3c1e --- /dev/null +++ b/stan/math/rev/fun/sylvester_solver.hpp @@ -0,0 +1,146 @@ +#ifndef STAN_MATH_REV_FUN_SYLVESTER_SOLVER_HPP +#define STAN_MATH_REV_FUN_SYLVESTER_SOLVER_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Solve the Sylvester equation A X + X B = C for X, reverse mode. + * + * Forward pass as in prim. With G the adjoint of X (real part: adjoint of + * Re X, imaginary part: adjoint of Im X) and Z the solution of + * A^H Z + Z B^H = G, the gradients are + * grad_A = -Z X^H, grad_B = -X^H Z, grad_C = Z. + * Complex inputs receive Re and Im of their gradient; real inputs receive + * only the real part. + * + * @tparam MatA type of A (real or complex, var or double scalars) + * @tparam MatB type of B (real or complex, var or double scalars) + * @tparam MatC type of C (real or complex, var or double scalars) + * @param[in] A n x n matrix + * @param[in] B m x m matrix + * @param[in] C n x m matrix + * @return complex n x m solution X with std::complex entries + * @throw std::invalid_argument if A or B is not square, or if the + * dimensions of C do not match + */ +template * = nullptr, + require_any_t>, is_var>, + is_var>>* = nullptr> +inline auto sylvester_solver(const MatA& A, const MatB& B, const MatC& C) { + using cdouble = std::complex; + using ret_type + = Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic>; + + check_square("sylvester_solver", "A", A); + check_square("sylvester_solver", "B", B); + check_size_match("sylvester_solver", "rows of A", A.rows(), "rows of C", + C.rows()); + check_size_match("sylvester_solver", "rows of B", B.rows(), "cols of C", + C.cols()); + if (unlikely(A.rows() == 0 || B.rows() == 0)) { + return ret_type(A.rows(), B.rows()); + } + + auto arena_A = to_arena(A); + auto arena_B = to_arena(B); + auto arena_C = to_arena(C); + + const Eigen::Index m = B.rows(); + + // ---- forward pass: identical to prim, on values ---- + Eigen::ComplexSchur> schur_A( + value_of(arena_A)); + auto arena_U = to_arena(schur_A.matrixU()); + auto arena_S = to_arena(schur_A.matrixT()); + + Eigen::ComplexSchur> schur_B( + value_of(arena_B)); + auto arena_V = to_arena(schur_B.matrixU()); + auto arena_T = to_arena(schur_B.matrixT()); + + // Y holds D = U^H C V initially, overwritten column by column + auto Y = (arena_U.adjoint() * value_of(arena_C) * arena_V).eval(); + // deep copy: copying an arena_matrix would alias arena_S + auto Sshift = arena_S.eval(); + for (Eigen::Index i = 0; i < m; ++i) { + Y.col(i).noalias() -= Y.leftCols(i) * arena_T.col(i).head(i); + Sshift.diagonal() = arena_S.diagonal().array() + arena_T(i, i); + Sshift.template triangularView().solveInPlace(Y.col(i)); + } + arena_t arena_X = arena_U * Y * arena_V.adjoint(); + + auto res = to_arena(arena_X.unaryExpr( + [](const cdouble& x) { return std::complex(x.real(), x.imag()); })); + + // ---- reverse pass ---- + reverse_pass_callback([arena_A, arena_B, arena_C, arena_U, arena_S, arena_V, + arena_T, arena_X, res]() mutable { + const Eigen::Index rows = arena_S.rows(); + const Eigen::Index cols = arena_T.rows(); + + // G = adj(Re X) + i * adj(Im X) + Eigen::MatrixXcd G(rows, cols); + for (Eigen::Index c = 0; c < cols; ++c) { + for (Eigen::Index r = 0; r < rows; ++r) { + G(r, c) = cdouble(res(r, c).real().adj(), res(r, c).imag().adj()); + } + } + + // W holds E = U^H G V initially; solve S^H W + W T^H = E, + // overwritten column by column from last to first + auto W = (arena_U.adjoint() * G * arena_V).eval(); + auto Sshift_r = arena_S.eval(); + for (Eigen::Index j = cols - 1; j >= 0; --j) { + const Eigen::Index k = cols - 1 - j; // already-solved columns + W.col(j).noalias() -= W.rightCols(k) * arena_T.row(j).tail(k).adjoint(); + Sshift_r.diagonal() = arena_S.diagonal().array() + arena_T(j, j); + Sshift_r.template triangularView().adjoint().solveInPlace( + W.col(j)); + } + const auto Z = (arena_U * W * arena_V.adjoint()).eval(); + + // Re/Im of the gradient go to the adjoints of Re/Im of the input; + // real inputs only receive Re(grad) + auto accumulate = [](auto& in, const Eigen::MatrixXcd& g) { + using in_t = std::decay_t; + if constexpr (is_complex>::value) { + for (Eigen::Index c = 0; c < g.cols(); ++c) { + for (Eigen::Index r = 0; r < g.rows(); ++r) { + in(r, c).real().adj() += g(r, c).real(); + in(r, c).imag().adj() += g(r, c).imag(); + } + } + } else { + in.adj() += g.real(); + } + }; + + if constexpr (is_var>::value) { + accumulate(arena_A, -(Z * arena_X.adjoint())); + } + if constexpr (is_var>::value) { + accumulate(arena_B, -(arena_X.adjoint() * Z)); + } + if constexpr (is_var>::value) { + accumulate(arena_C, Z); + } + }); + // copy out of the arena: the callback reads adjoints from res, so the + // caller must not be able to overwrite its entries + return ret_type(res); +} + +} // namespace math +} // namespace stan +#endif \ No newline at end of file diff --git a/test/unit/math/mix/fun/sylvester_test.cpp b/test/unit/math/mix/fun/sylvester_test.cpp new file mode 100644 index 00000000000..7db8f7b7c00 --- /dev/null +++ b/test/unit/math/mix/fun/sylvester_test.cpp @@ -0,0 +1,147 @@ +#include + +namespace sylvester_solver_test { + +// Real scalar depending on both Re X and Im X, so both carry adjoints +template +stan::math::var weighted_sum(const MatX& X, const Eigen::MatrixXd& W_re, + const Eigen::MatrixXd& W_im) { + stan::math::var result = 0; + for (Eigen::Index c = 0; c < X.cols(); ++c) { + for (Eigen::Index r = 0; r < X.rows(); ++r) { + result += W_re(r, c) * X(r, c).real() + W_im(r, c) * X(r, c).imag(); + } + } + return result; +} + +} // namespace sylvester_solver_test + +TEST(mathMixMatFun, sylvester_solver) { + auto f = [](const auto& a, const auto& b, const auto& c) { + return stan::math::sylvester_solver(a, b, c); + }; + + // A has a complex eigenvalue pair, so the complex Schur path is used + Eigen::MatrixXd A(2, 2); + A << 4, 1, -1, 3; + Eigen::MatrixXd B(3, 3); + B << 5, 0, 1, 2, 6, 0, -1, 1, 4; + Eigen::MatrixXd C(2, 3); + C << 1, 2, 3, 4, 5, 6; + stan::test::expect_ad(f, A, B, C); + // expect_ad_matvar does not apply: the result has complex scalars, + // which have no var_value form. See sylvester_solver_varmat. + + // Lyapunov: B = A^T shares A's autodiff variables + auto g = [](const auto& a, const auto& c) { + return stan::math::sylvester_solver(a, a.transpose(), c); + }; + Eigen::MatrixXd Q(2, 2); + Q << 1, 2, 3, 4; + stan::test::expect_ad(g, A, Q); + + // size mismatch must throw for every autodiff type, as it does for double + Eigen::MatrixXd C_bad(3, 3); + C_bad.setOnes(); + stan::test::expect_ad(f, A, B, C_bad); +} + +TEST(mathMixMatFun, sylvester_solver_complex) { + auto f = [](const auto& a, const auto& b, const auto& c) { + return stan::math::sylvester_solver(a, b, c); + }; + using cdouble = std::complex; + + // eigenvalue sums lambda_i(A) + mu_j(B) have real part around 9 + Eigen::MatrixXcd A(2, 2); + A << cdouble(4, 1), cdouble(1, -0.5), cdouble(-1, 0.3), cdouble(3, -0.2); + Eigen::MatrixXcd B(2, 2); + B << cdouble(5, -0.4), cdouble(0.5, 1), cdouble(2, -0.7), cdouble(6, 0.3); + Eigen::MatrixXcd C(2, 2); + C << cdouble(1, 2), cdouble(-0.5, 1), cdouble(3, -1), cdouble(0.7, 0.2); + stan::test::expect_ad(f, A, B, C); + + // real and complex arguments mixed + Eigen::MatrixXd A_r = A.real(); + Eigen::MatrixXd B_r = B.real(); + Eigen::MatrixXd C_r = C.real(); + // real system, complex right-hand side + stan::test::expect_ad(f, A_r, B_r, C); + // complex A only + stan::test::expect_ad(f, A, B_r, C_r); +} + +TEST(mathMixMatFun, sylvester_solver_varmat) { + // var_value inputs must give the same value and adjoints as + // Matrix inputs + using stan::math::var; + using stan::math::var_value; + using sylvester_solver_test::weighted_sum; + Eigen::MatrixXd A(2, 2); + A << 4, 1, -1, 3; + Eigen::MatrixXd B(3, 3); + B << 5, 0, 1, 2, 6, 0, -1, 1, 4; + Eigen::MatrixXd C(2, 3); + C << 1, 2, 3, 4, 5, 6; + Eigen::MatrixXd W_re(2, 3); + W_re << 0.3, -1.2, 0.7, 2.1, -0.4, 1.5; + Eigen::MatrixXd W_im(2, 3); + W_im << -0.8, 0.6, 1.9, -1.1, 0.2, -0.5; + + var_value A_vm(A); + var_value B_vm(B); + var_value C_vm(C); + var f_vm = weighted_sum(stan::math::sylvester_solver(A_vm, B_vm, C_vm), W_re, + W_im); + f_vm.grad(); + const double val_vm = f_vm.val(); + Eigen::MatrixXd adj_A_vm = A_vm.adj(); + Eigen::MatrixXd adj_B_vm = B_vm.adj(); + Eigen::MatrixXd adj_C_vm = C_vm.adj(); + stan::math::recover_memory(); + + auto A_mv = stan::math::to_var(A); + auto B_mv = stan::math::to_var(B); + auto C_mv = stan::math::to_var(C); + var f_mv = weighted_sum(stan::math::sylvester_solver(A_mv, B_mv, C_mv), W_re, + W_im); + f_mv.grad(); + + stan::test::expect_near_rel("value", val_vm, f_mv.val()); + stan::test::expect_near_rel("adjoint of A", adj_A_vm, + Eigen::MatrixXd(A_mv.adj())); + stan::test::expect_near_rel("adjoint of B", adj_B_vm, + Eigen::MatrixXd(B_mv.adj())); + stan::test::expect_near_rel("adjoint of C", adj_C_vm, + Eigen::MatrixXd(C_mv.adj())); + stan::math::recover_memory(); +} + +TEST(mathMixMatFun, sylvester_solver_result_is_independent) { + // Writing into the result must not change what the reverse pass reads: + // after X(0, 0) = 2 X(0, 0), d Re X(0, 0) / dC must double + using stan::math::var; + Eigen::MatrixXd A(2, 2); + A << 4, 1, -1, 3; + Eigen::MatrixXd B(2, 2); + B << 5, 0, 2, 6; + Eigen::MatrixXd C(2, 2); + C << 1, 2, 3, 4; + + auto C1 = stan::math::to_var(C); + auto X1 = stan::math::sylvester_solver(A, B, C1); + var f1 = X1(0, 0).real(); + f1.grad(); + Eigen::MatrixXd adj_C = C1.adj(); + stan::math::recover_memory(); + + auto C2 = stan::math::to_var(C); + auto X2 = stan::math::sylvester_solver(A, B, C2); + X2(0, 0) = 2.0 * X2(0, 0); + var f2 = X2(0, 0).real(); + f2.grad(); + stan::test::expect_near_rel("adjoint of C", Eigen::MatrixXd(2.0 * adj_C), + Eigen::MatrixXd(C2.adj())); + stan::math::recover_memory(); +} From 227b4413522d19ec6862c4c7bd6e5b6fa0d9c8bf Mon Sep 17 00:00:00 2001 From: Stefan Kienle <127690239+stekien@users.noreply.github.com> Date: Thu, 15 Oct 2026 17:51:39 -0400 Subject: [PATCH 4/6] add zero-size prim test for sylvester solver Co-Authored-By: Claude Opus 5.5 --- test/unit/math/prim/fun/sylvester_test.cpp | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/test/unit/math/prim/fun/sylvester_test.cpp b/test/unit/math/prim/fun/sylvester_test.cpp index 2e49bd24386..5f77c2d8590 100644 --- a/test/unit/math/prim/fun/sylvester_test.cpp +++ b/test/unit/math/prim/fun/sylvester_test.cpp @@ -155,4 +155,18 @@ TEST(MathFunctions, sylvester_solver_errors) { // cols of C do not match B Eigen::MatrixXd C_cols = Eigen::MatrixXd::Ones(3, 3); EXPECT_THROW(sylvester_solver(A, B, C_cols), std::invalid_argument); -} \ No newline at end of file +} + +TEST(MathFunctions, sylvester_solver_zero_size) { + using stan::math::sylvester_solver; + Eigen::MatrixXd I2 = Eigen::MatrixXd::Identity(2, 2); + Eigen::MatrixXd E(0, 0); + + Eigen::MatrixXcd X1 = sylvester_solver(E, I2, Eigen::MatrixXd(0, 2)); + EXPECT_EQ(X1.rows(), 0); + EXPECT_EQ(X1.cols(), 2); + + Eigen::MatrixXcd X2 = sylvester_solver(I2, E, Eigen::MatrixXd(2, 0)); + EXPECT_EQ(X2.rows(), 2); + EXPECT_EQ(X2.cols(), 0); +} From 42ce8c6856780cc9e435b3e3bfb943798d78edc6 Mon Sep 17 00:00:00 2001 From: Stefan Kienle <127690239+stekien@users.noreply.github.com> Date: Fri, 16 Oct 2026 09:32:55 -0400 Subject: [PATCH 5/6] fix formatting --- stan/math/prim/fun/sylvester_solver.hpp | 38 ++++++++++------------ stan/math/rev/fun/sylvester_solver.hpp | 4 +-- test/unit/math/prim/fun/sylvester_test.cpp | 12 +++---- 3 files changed, 26 insertions(+), 28 deletions(-) diff --git a/stan/math/prim/fun/sylvester_solver.hpp b/stan/math/prim/fun/sylvester_solver.hpp index ebac314509a..4ee68ade2ed 100644 --- a/stan/math/prim/fun/sylvester_solver.hpp +++ b/stan/math/prim/fun/sylvester_solver.hpp @@ -29,29 +29,27 @@ namespace math { * @throw std::invalid_argument if A or B is not square, or if the * dimensions of C do not match */ -template * = nullptr, - require_all_not_t>, - is_var>, - is_var>>* = nullptr> -inline auto sylvester_solver(const MatA& A, const MatB& B, - const MatC& C) { - +template < + typename MatA, typename MatB, typename MatC, + require_all_matrix_t* = nullptr, + require_all_not_t>, is_var>, + is_var>>* = nullptr> +inline auto sylvester_solver(const MatA& A, const MatB& B, const MatC& C) { // common base scalar: double, or fvar<...> if any input is fvar - using base_t = return_type_t, base_type_t, - base_type_t>; + using base_t + = return_type_t, base_type_t, base_type_t>; using complex_t = std::complex; using complex_mat_t = Eigen::Matrix; // A and B are promoted to base_t but stay real or complex as given - using A_mat_t = Eigen::Matrix< - std::conditional_t>::value, complex_t, - base_t>, - Eigen::Dynamic, Eigen::Dynamic>; - using B_mat_t = Eigen::Matrix< - std::conditional_t>::value, complex_t, - base_t>, - Eigen::Dynamic, Eigen::Dynamic>; + using A_mat_t + = Eigen::Matrix>::value, + complex_t, base_t>, + Eigen::Dynamic, Eigen::Dynamic>; + using B_mat_t + = Eigen::Matrix>::value, + complex_t, base_t>, + Eigen::Dynamic, Eigen::Dynamic>; check_square("sylvester_solver", "A", A); check_square("sylvester_solver", "B", B); @@ -60,7 +58,7 @@ inline auto sylvester_solver(const MatA& A, const MatB& B, check_size_match("sylvester_solver", "rows of B", B.rows(), "cols of C", C.cols()); - if(unlikely(A.rows() == 0 || B.rows() == 0)) { + if (unlikely(A.rows() == 0 || B.rows() == 0)) { return complex_mat_t(A.rows(), B.rows()); } @@ -93,4 +91,4 @@ inline auto sylvester_solver(const MatA& A, const MatB& B, } // namespace math } // namespace stan -#endif \ No newline at end of file +#endif diff --git a/stan/math/rev/fun/sylvester_solver.hpp b/stan/math/rev/fun/sylvester_solver.hpp index 539e7aa3c1e..ee1b8175123 100644 --- a/stan/math/rev/fun/sylvester_solver.hpp +++ b/stan/math/rev/fun/sylvester_solver.hpp @@ -106,7 +106,7 @@ inline auto sylvester_solver(const MatA& A, const MatB& B, const MatC& C) { W.col(j).noalias() -= W.rightCols(k) * arena_T.row(j).tail(k).adjoint(); Sshift_r.diagonal() = arena_S.diagonal().array() + arena_T(j, j); Sshift_r.template triangularView().adjoint().solveInPlace( - W.col(j)); + W.col(j)); } const auto Z = (arena_U * W * arena_V.adjoint()).eval(); @@ -143,4 +143,4 @@ inline auto sylvester_solver(const MatA& A, const MatB& B, const MatC& C) { } // namespace math } // namespace stan -#endif \ No newline at end of file +#endif diff --git a/test/unit/math/prim/fun/sylvester_test.cpp b/test/unit/math/prim/fun/sylvester_test.cpp index 5f77c2d8590..5fc5c4b5f62 100644 --- a/test/unit/math/prim/fun/sylvester_test.cpp +++ b/test/unit/math/prim/fun/sylvester_test.cpp @@ -110,18 +110,18 @@ TEST(MathFunctions, sylvester_solver_lyapunov) { EXPECT_LT((Xr - Xr.transpose()).norm(), 1e-12); // stable A and positive definite Q give a positive definite solution - Eigen::SelfAdjointEigenSolver eig( - 0.5 * (Xr + Xr.transpose())); + Eigen::SelfAdjointEigenSolver eig(0.5 + * (Xr + Xr.transpose())); EXPECT_GT(eig.eigenvalues().minCoeff(), 0.0); } TEST(MathFunctions, sylvester_solver_expression_inputs) { // Eigen expressions go through plain_type_t; result must match plain inputs using stan::math::sylvester_solver; - Eigen::MatrixXd A = Eigen::MatrixXd::Random(4, 4) - + 8.0 * Eigen::MatrixXd::Identity(4, 4); - Eigen::MatrixXd B = Eigen::MatrixXd::Random(3, 3) - + 6.0 * Eigen::MatrixXd::Identity(3, 3); + Eigen::MatrixXd A + = Eigen::MatrixXd::Random(4, 4) + 8.0 * Eigen::MatrixXd::Identity(4, 4); + Eigen::MatrixXd B + = Eigen::MatrixXd::Random(3, 3) + 6.0 * Eigen::MatrixXd::Identity(3, 3); Eigen::MatrixXd C = Eigen::MatrixXd::Random(4, 3); Eigen::MatrixXd At = A.transpose(); From ddad81072c5f0f22ff2a6b0fd84587e3fed4e5cb Mon Sep 17 00:00:00 2001 From: Stefan Kienle <127690239+stekien@users.noreply.github.com> Date: Fri, 16 Oct 2026 16:19:20 -0400 Subject: [PATCH 6/6] fix compiler error ofr fvar type in prim/sylvester --- stan/math/prim/fun/sylvester_solver.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/prim/fun/sylvester_solver.hpp b/stan/math/prim/fun/sylvester_solver.hpp index 4ee68ade2ed..bb59114dc9a 100644 --- a/stan/math/prim/fun/sylvester_solver.hpp +++ b/stan/math/prim/fun/sylvester_solver.hpp @@ -76,7 +76,8 @@ inline auto sylvester_solver(const MatA& A, const MatB& B, const MatC& C) { const auto& T = schur_B.matrixT(); // Y holds D = U^H C V initially, overwritten column by column - auto Y = (U.adjoint() * C * V).eval(); + const auto C_c = C.template cast(); + auto Y = (U.adjoint() * C_c * V).eval(); auto Sshift = S; for (Eigen::Index i = 0; i < m; ++i) { Y.col(i).noalias() -= Y.leftCols(i) * T.col(i).head(i);