From d08eaad37b7f5808bd423f951773167ce2abee79 Mon Sep 17 00:00:00 2001 From: Stefan Kienle <127690239+stekien@users.noreply.github.com> Date: Tue, 3 Nov 2026 13:07:49 -0500 Subject: [PATCH] add solve_lyapunov function as a special case of the solve_sylvester with respective tests --- stan/math/fwd/fun/solve_sylvester.hpp | 81 +++++++++++ stan/math/prim/fun/solve_sylvester.hpp | 56 ++++++++ stan/math/rev/fun/solve_sylvester.hpp | 117 ++++++++++++++++ .../math/mix/fun/solve_sylvester_test.cpp | 62 +++++++-- .../math/prim/fun/solve_sylvester_test.cpp | 128 ++++++++++++++---- 5 files changed, 410 insertions(+), 34 deletions(-) diff --git a/stan/math/fwd/fun/solve_sylvester.hpp b/stan/math/fwd/fun/solve_sylvester.hpp index 644cb26270a..58f1fc5fc93 100644 --- a/stan/math/fwd/fun/solve_sylvester.hpp +++ b/stan/math/fwd/fun/solve_sylvester.hpp @@ -103,6 +103,87 @@ inline auto solve_sylvester(const MatA& A, const MatB& B, const MatC& C) { return res; } +/** + * Solve the continuous Lyapunov equation A X + X A^H = C for X, forward mode. + * + * The value X solves the equation on the values of A and C. Differentiating + * the equation gives the tangent as the solution of another Lyapunov equation + * with the same coefficient, + * A X' + X' A^H = C' - A' X - X A'^H. + * Both solves run on the inner scalar type, so no Schur decomposition is taken + * on fvar values. + * + * @tparam MatA type of A (real or complex, fvar or arithmetic scalars) + * @tparam MatC type of C (real or complex, fvar or arithmetic scalars) + * @param[in] A n x n matrix + * @param[in] C n x n matrix + * @return complex n x n solution X with std::complex entries + * @throw std::invalid_argument if A is not square, or if the dimensions of C + * do not match + */ +template * = nullptr, + require_any_bt_fvar* = nullptr> +inline auto solve_lyapunov(const MatA& A, const MatC& C) { + using fvar_t = return_type_t, base_type_t>; + using ret_type + = Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic>; + + check_square("solve_lyapunov", "A", A); + check_size_match("solve_lyapunov", "rows of A", A.rows(), "rows of C", + C.rows()); + check_size_match("solve_lyapunov", "rows of A", A.rows(), "cols of C", + C.cols()); + if (unlikely(A.rows() == 0)) { + return ret_type(0, 0); + } + + using complex_inner_t = std::complex; + + // tangent of an fvar or complex matrix, staying real or complex + auto tangent = [](const auto& x) { + return x + .unaryExpr([](const auto& v) { + if constexpr (is_complex>::value) { + return std::complex>( + v.real().d_, v.imag().d_); + } else { + return v.d_; + } + }) + .eval(); + }; + + const auto& A_ref = to_ref(A); + const auto& C_ref = to_ref(C); + const auto& A_val = to_ref(value_of(A_ref)); + using inner_mat_t + = Eigen::Matrix; + + auto X = solve_lyapunov(A_val, value_of(C_ref)); + + inner_mat_t rhs = inner_mat_t::Zero(X.rows(), X.cols()); + if constexpr (is_fvar>::value) { + rhs += tangent(C_ref); + } + if constexpr (is_fvar>::value) { + const auto A_d = tangent(A_ref); + rhs -= A_d * X + X * A_d.adjoint(); + } + inner_mat_t X_d = solve_lyapunov(A_val, rhs); + const Eigen::Index X_rows = X.rows(); + const Eigen::Index X_cols = X.cols(); + ret_type res(X_rows, X_cols); + for (Eigen::Index j = 0; j < X_cols; ++j) { + for (Eigen::Index i = 0; i < X_rows; ++i) { + res(i, j) + = std::complex(fvar_t(X(i, j).real(), X_d(i, j).real()), + fvar_t(X(i, j).imag(), X_d(i, j).imag())); + } + } + return res; +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/prim/fun/solve_sylvester.hpp b/stan/math/prim/fun/solve_sylvester.hpp index b1b14d6c2a1..544f94046e5 100644 --- a/stan/math/prim/fun/solve_sylvester.hpp +++ b/stan/math/prim/fun/solve_sylvester.hpp @@ -70,6 +70,62 @@ inline auto solve_sylvester(MatA&& A, MatB&& B, MatC&& C) { return complex_mat_t(U * Y * V.adjoint()); } +/** + * Solve the continuous Lyapunov equation A X + X A^H = C for X. + * + * Special case of solve_sylvester with B = A^H. The complex Schur + * decomposition A = U T U^H also gives A^H = U T^H U^H, so only one + * decomposition is needed. The equation reduces to T Y + Y T^H = U^H C U, + * solved column by column from last to first, and X = U Y U^H. A unique + * solution requires lambda_i(A) + conj(lambda_j(A)) != 0 for all + * eigenvalues of A. + * + * @tparam MatA type of A (real or complex arithmetic scalars) + * @tparam MatC type of C (real or complex arithmetic scalars) + * @param[in] A n x n matrix + * @param[in] C n x n matrix + * @return complex n x n solution X + * @throw std::invalid_argument if A is not square, or if the dimensions of C + * do not match + */ +template * = nullptr, + require_all_bt_arithmetic* = nullptr> +inline auto solve_lyapunov(MatA&& A, MatC&& C) { + using complex_mat_t + = Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic>; + + check_square("solve_lyapunov", "A", A); + check_size_match("solve_lyapunov", "rows of A", A.rows(), "rows of C", + C.rows()); + check_size_match("solve_lyapunov", "rows of A", A.rows(), "cols of C", + C.cols()); + + if (unlikely(A.rows() == 0)) { + return complex_mat_t(0, 0); + } + + const Eigen::Index n = A.rows(); + + // ComplexSchur needs a plain matrix (same pattern as complex_schur_decompose) + Eigen::ComplexSchur> schur_A{eval(std::forward(A))}; + auto&& U = schur_A.matrixU(); + auto&& T = schur_A.matrixT(); + + // Y holds D = U^H C U initially; T^H is lower triangular, so Y is + // overwritten column by column from last to first + auto&& C_ref = to_ref(std::forward(C)); + auto Y = (U.adjoint() * C_ref * U).eval(); + plain_type_t Tshift = T; + for (Eigen::Index j = n - 1; j >= 0; --j) { + const Eigen::Index k = n - 1 - j; // already-solved columns + Y.col(j).noalias() -= Y.rightCols(k) * T.row(j).tail(k).adjoint(); + Tshift.diagonal() = T.diagonal().array() + std::conj(T(j, j)); + Tshift.template triangularView().solveInPlace(Y.col(j)); + } + return complex_mat_t(U * Y * U.adjoint()); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/solve_sylvester.hpp b/stan/math/rev/fun/solve_sylvester.hpp index 0e6c968b508..f4592aa3b6f 100644 --- a/stan/math/rev/fun/solve_sylvester.hpp +++ b/stan/math/rev/fun/solve_sylvester.hpp @@ -138,6 +138,123 @@ inline auto solve_sylvester(MatA&& A, MatB&& B, MatC&& C) { return ret_type(res); } +/** + * Solve the continuous Lyapunov equation A X + X A^H = C for X, reverse mode. + * + * Forward pass as in prim, with a single complex Schur decomposition + * A = U T U^H. 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 A = G, + * the gradients are + * grad_A = -(Z X^H + Z^H X), grad_C = Z, + * where the second term of grad_A comes from A^H in place of B. The adjoint + * equation reuses the same Schur factors as T^H W + W T = U^H G U. + * 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 MatC type of C (real or complex, var or double scalars) + * @param[in] A n x n matrix + * @param[in] C n x n matrix + * @return complex n x n solution X with std::complex entries + * @throw std::invalid_argument if A is not square, or if the dimensions of C + * do not match + */ +template * = nullptr, + require_any_bt_var* = nullptr> +inline auto solve_lyapunov(MatA&& A, MatC&& C) { + using cdouble = std::complex; + using ret_type + = Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic>; + + check_square("solve_lyapunov", "A", A); + check_size_match("solve_lyapunov", "rows of A", A.rows(), "rows of C", + C.rows()); + check_size_match("solve_lyapunov", "rows of A", A.rows(), "cols of C", + C.cols()); + if (unlikely(A.rows() == 0)) { + return ret_type(0, 0); + } + + auto arena_A = to_arena(std::forward(A)); + auto arena_C = to_arena(std::forward(C)); + + const Eigen::Index n = A.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_T = to_arena(schur_A.matrixT()); + + // Y holds D = U^H C U initially, overwritten column by column from last + // to first + auto Y = (arena_U.adjoint() * value_of(arena_C) * arena_U).eval(); + plain_type_t Tshift = arena_T; + for (Eigen::Index j = n - 1; j >= 0; --j) { + const Eigen::Index k = n - 1 - j; // already-solved columns + Y.col(j).noalias() -= Y.rightCols(k) * arena_T.row(j).tail(k).adjoint(); + Tshift.diagonal() = arena_T.diagonal().array() + std::conj(arena_T(j, j)); + Tshift.template triangularView().solveInPlace(Y.col(j)); + } + arena_t arena_X = arena_U * Y * arena_U.adjoint(); + + auto res = to_arena(arena_X.cast>()); + + // ---- reverse pass ---- + reverse_pass_callback([arena_A, arena_C, arena_U, arena_T, arena_X, + res]() mutable { + const Eigen::Index n = arena_T.rows(); + + // G = adj(Re X) + i * adj(Im X) + Eigen::MatrixXcd G(n, n); + for (Eigen::Index c = 0; c < n; ++c) { + for (Eigen::Index r = 0; r < n; ++r) { + G(r, c) = cdouble(res(r, c).real().adj(), res(r, c).imag().adj()); + } + } + + // W holds E = U^H G U initially; solve T^H W + W T = E, + // overwritten column by column from first to last + auto W = (arena_U.adjoint() * G * arena_U).eval(); + auto Tshift_r = arena_T.eval(); + for (Eigen::Index j = 0; j < n; ++j) { + W.col(j).noalias() -= W.leftCols(j) * arena_T.col(j).head(j); + Tshift_r.diagonal() + = arena_T.diagonal().array() + std::conj(arena_T(j, j)); + Tshift_r.template triangularView().adjoint().solveInPlace( + W.col(j)); + } + const auto Z = (arena_U * W * arena_U.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) { + const Eigen::Index rows = g.rows(); + const Eigen::Index cols = g.cols(); + for (Eigen::Index c = 0; c < cols; ++c) { + for (Eigen::Index r = 0; r < 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() + Z.adjoint() * arena_X)); + } + if constexpr (is_var>::value) { + accumulate(arena_C, Z); + } + }); + return ret_type(res); +} + } // namespace math } // namespace stan #endif diff --git a/test/unit/math/mix/fun/solve_sylvester_test.cpp b/test/unit/math/mix/fun/solve_sylvester_test.cpp index 7683220adce..f7e52d7352f 100644 --- a/test/unit/math/mix/fun/solve_sylvester_test.cpp +++ b/test/unit/math/mix/fun/solve_sylvester_test.cpp @@ -33,14 +33,6 @@ TEST(mathMixMatFun, solve_sylvester) { // expect_ad_matvar does not apply: the result has complex scalars, // which have no var_value form. See solve_sylvester_varmat. - // Lyapunov: B = A^T shares A's autodiff variables - auto g = [](const auto& a, const auto& c) { - return stan::math::solve_sylvester(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(); @@ -164,3 +156,57 @@ TEST(mathMixMatFun, solve_sylvester_result_is_independent) { Eigen::MatrixXd(C2.adj())); stan::math::recover_memory(); } + +TEST(mathMixMatFun, solve_lyapunov) { + auto f = [](const auto& a, const auto& c) { + return stan::math::solve_lyapunov(a, c); + }; + + // A has a complex eigenvalue pair, so the complex Schur path is used + Eigen::MatrixXd A(3, 3); + A << -1, -1, 0, 1, -1, 0.5, 0.3, -0.2, -2; + Eigen::MatrixXd C(3, 3); + C << 1, 2, 3, 4, 5, 6, 7, 8, 9; + stan::test::expect_ad(f, A, C); + + // size mismatch must throw for every autodiff type, as it does for double + Eigen::MatrixXd C_bad(3, 2); + C_bad.setOnes(); + EXPECT_THROW(stan::math::solve_lyapunov(A, C_bad), std::invalid_argument); + + auto A_v = stan::math::to_var(A); + auto C_bad_v = stan::math::to_var(C_bad); + EXPECT_THROW(stan::math::solve_lyapunov(A_v, C_bad_v), std::invalid_argument); + stan::math::recover_memory(); +} + +TEST(mathMixMatFun, solve_lyapunov_complex) { + auto f = [](const auto& a, const auto& c) { + return stan::math::solve_lyapunov(a, c); + }; + using cdouble = std::complex; + + // eigenvalue sums lambda_i + conj(lambda_j) have real part around -7 + Eigen::MatrixXcd A(2, 2); + A << cdouble(-4, 1), cdouble(1, -0.5), cdouble(-1, 0.3), cdouble(-3, -0.2); + 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, C); + + // real system, complex right-hand side + Eigen::MatrixXd A_r = A.real(); + stan::test::expect_ad(f, A_r, C); +} + +TEST(mathMixMatFun, solve_lyapunov_triangular) { + // Schur of a triangular matrix deflates on exact zeros, which must not + // drop the tangents of the zero entries + auto f = [](const auto& a, const auto& c) { + return stan::math::solve_lyapunov(a, c); + }; + Eigen::MatrixXd A(3, 3); + A << -5, 2, -1, 0, -6, 1, 0, 0, -4; + Eigen::MatrixXd C(3, 3); + C << 1, 2, 3, 4, 5, 6, 7, 8, 9; + stan::test::expect_ad(f, A, C); +} diff --git a/test/unit/math/prim/fun/solve_sylvester_test.cpp b/test/unit/math/prim/fun/solve_sylvester_test.cpp index 2d70cfe1b90..d449c274009 100644 --- a/test/unit/math/prim/fun/solve_sylvester_test.cpp +++ b/test/unit/math/prim/fun/solve_sylvester_test.cpp @@ -89,32 +89,6 @@ TEST(MathFunctions, solve_sylvester_diagonal) { } } -TEST(MathFunctions, solve_sylvester_lyapunov) { - // continuous Lyapunov equation A X + X A^T = -Q - using stan::math::solve_sylvester; - 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 = solve_sylvester(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, solve_sylvester_expression_inputs) { // Eigen expressions go through plain_type_t; result must match plain inputs using stan::math::solve_sylvester; @@ -170,3 +144,105 @@ TEST(MathFunctions, solve_sylvester_zero_size) { EXPECT_EQ(X2.rows(), 2); EXPECT_EQ(X2.cols(), 0); } + +// Stable A (by Gershgorin all eigenvalues have Re <= -(n + 1) < 0) keeps +// lambda_i + conj(lambda_j) away from zero. + +TEST(MathFunctions, solve_lyapunov_real_known_solution) { + using stan::math::solve_lyapunov; + const int n = 6; + Eigen::MatrixXd A = Eigen::MatrixXd::Random(n, n) + - 2.0 * n * Eigen::MatrixXd::Identity(n, n); + Eigen::MatrixXd X_true = Eigen::MatrixXd::Random(n, n); + Eigen::MatrixXd C = A * X_true + X_true * A.transpose(); + + Eigen::MatrixXcd X = solve_lyapunov(A, C); + + EXPECT_EQ(X.rows(), n); + EXPECT_EQ(X.cols(), n); + // 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 * A.transpose() - C).norm(), 1e-10); +} + +TEST(MathFunctions, solve_lyapunov_complex_known_solution) { + using stan::math::solve_lyapunov; + const int n = 6; + Eigen::MatrixXcd A = Eigen::MatrixXcd::Random(n, n) + - 2.0 * n * Eigen::MatrixXcd::Identity(n, n); + Eigen::MatrixXcd X_true = Eigen::MatrixXcd::Random(n, n); + Eigen::MatrixXcd C = A * X_true + X_true * A.adjoint(); + + Eigen::MatrixXcd X = solve_lyapunov(A, C); + + EXPECT_LT((X - X_true).norm(), 1e-10); + EXPECT_LT((A * X + X * A.adjoint() - C).norm(), 1e-10); +} + +TEST(MathFunctions, solve_lyapunov_matches_sylvester) { + using stan::math::solve_lyapunov; + using stan::math::solve_sylvester; + const int n = 5; + Eigen::MatrixXcd A = Eigen::MatrixXcd::Random(n, n) + - 2.0 * n * Eigen::MatrixXcd::Identity(n, n); + Eigen::MatrixXcd C = Eigen::MatrixXcd::Random(n, n); + Eigen::MatrixXcd A_adj = A.adjoint(); + + Eigen::MatrixXcd X_lyap = solve_lyapunov(A, C); + Eigen::MatrixXcd X_sylv = solve_sylvester(A, A_adj, C); + + EXPECT_LT((X_lyap - X_sylv).norm(), 1e-12); + + // real A, complex C + Eigen::MatrixXd A_r = A.real(); + Eigen::MatrixXd A_rt = A_r.transpose(); + EXPECT_LT((solve_lyapunov(A_r, C) - solve_sylvester(A_r, A_rt, C)).norm(), + 1e-12); +} + +TEST(MathFunctions, solve_lyapunov_hermitian) { + // stable A and Hermitian negative definite C = -Q give a Hermitian + // positive definite solution + using stan::math::solve_lyapunov; + const int n = 5; + Eigen::MatrixXcd A = Eigen::MatrixXcd::Random(n, n) + - 2.0 * n * Eigen::MatrixXcd::Identity(n, n); + Eigen::MatrixXcd L = Eigen::MatrixXcd::Random(n, n); + Eigen::MatrixXcd Q = L * L.adjoint() + Eigen::MatrixXcd::Identity(n, n); + + Eigen::MatrixXcd X = solve_lyapunov(A, -Q); + + EXPECT_LT((A * X + X * A.adjoint() + Q).norm(), 1e-10); + EXPECT_LT((X - X.adjoint()).norm(), 1e-12); + Eigen::SelfAdjointEigenSolver eig(0.5 * (X + X.adjoint())); + EXPECT_GT(eig.eigenvalues().minCoeff(), 0.0); +} + +TEST(MathFunctions, solve_lyapunov_errors) { + using stan::math::solve_lyapunov; + Eigen::MatrixXd A = Eigen::MatrixXd::Identity(3, 3); + Eigen::MatrixXd C = Eigen::MatrixXd::Ones(3, 3); + + // valid input does not throw + EXPECT_NO_THROW(solve_lyapunov(A, C)); + + // A not square + Eigen::MatrixXd A_rect = Eigen::MatrixXd::Ones(3, 2); + EXPECT_THROW(solve_lyapunov(A_rect, C), std::invalid_argument); + + // rows of C do not match A + Eigen::MatrixXd C_rows = Eigen::MatrixXd::Ones(4, 3); + EXPECT_THROW(solve_lyapunov(A, C_rows), std::invalid_argument); + + // cols of C do not match A + Eigen::MatrixXd C_cols = Eigen::MatrixXd::Ones(3, 4); + EXPECT_THROW(solve_lyapunov(A, C_cols), std::invalid_argument); +} + +TEST(MathFunctions, solve_lyapunov_zero_size) { + Eigen::MatrixXd E(0, 0); + Eigen::MatrixXcd X = stan::math::solve_lyapunov(E, E); + EXPECT_EQ(X.rows(), 0); + EXPECT_EQ(X.cols(), 0); +}