[stekien]

Summary

Adds solve_lyapunov(A, C), which solves the continuous Lyapunov equation A X + X A^H = C. For real A this is the usual A X + X A^T = C. It is a special case of solve_sylvester(A, B, C) (added in #3467) with B = A^H. It has its own implementation so it needs only one Schur decomposition.

Depends on #3499. This branch is built on the rename of sylvester_solver to solve_sylvester in #3499 and includes that commit. Please merge #3499 first. After that, the diff here is only the Lyapunov commit.

The code is in stan/math/{prim,rev,fwd}/fun/solve_sylvester.hpp, next to solve_sylvester:

• prim: Computes one complex Schur decomposition A = U T U^H, which also gives A^H = U T^H U^H. The equation becomes T Y + Y T^H = U^H C U. This is solved one column at a time, from last to first, with a triangular solve per column. Then X = U Y U^H.
• rev: The forward pass is the same as prim, on values. U, T and X are stored in the arena. In the reverse pass, with G the adjoint of X, it solves the adjoint equation A^H Z + Z A = G with the same Schur factors. It then adds grad_A = -(Z X^H + Z^H X) and grad_C = Z. Complex inputs receive the real and imaginary parts of the gradient; real inputs receive only the real part.
• fwd: Differentiating the equation gives A X' + X' A^H = C' - A' X - X A'^H, which is another Lyapunov equation with the same A. Both solves run on the inner scalar type, so no Schur decomposition is taken on fvar values.

Behavior shared with solve_sylvester:

• Inputs can be real or complex.
• The result is always a complex matrix.
• A unique solution requires lambda_i(A) + conj(lambda_j(A)) != 0 for all pairs of eigenvalues of A.
• It throws std::invalid_argument if A is not square or if C is not the same size as A.

Tests

prim (test/unit/math/prim/fun/solve_sylvester_test.cpp):

• Real and complex inputs with known solutions. Checks the solution error and the residual. Real inputs give a real result up to roundoff.
• The result matches solve_sylvester(A, A^H, C), both for complex A and for real A with complex C.
• A stable A and a Hermitian negative definite C give a Hermitian positive definite solution.
• Errors for a non-square A and for C with the wrong number of rows or columns.
• Zero-size input.

mix (test/unit/math/mix/fun/solve_sylvester_test.cpp), using expect_ad:

• A real 3x3 A with a pair of complex eigenvalues.
• Complex A and C, and real A with complex C.
• A triangular A. Here the Schur decomposition deflates on exact zeros, and the test checks that the tangents of those zero entries are kept.
• A size mismatch throws for both double and var inputs.

These replace the earlier Lyapunov tests that called solve_sylvester(A, A.transpose(), C).

Side Effects

None, apart from the merge order with #3499.

Release notes

Added solve_lyapunov(A, C), which solves the continuous Lyapunov equation AX + XA^H = C for real or complex square matrices, with forward- and reverse-mode autodiff.

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

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.25	0.24	1.02	1.85% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	0.99	-1.1% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.96	-4.03% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.69	0.69	1.0	0.39% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.41	3.4	1.0	0.14% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	42.9	42.14	1.02	1.76% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.6	0.6	1.0	-0.35% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.66	4.75	0.98	-1.97% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	22.23	21.73	1.02	2.23% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	1.01	1.01	0.99	-0.69% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	0.44% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	9.32	9.33	1.0	-0.11% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.59	6.62	1.0	-0.5% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	229.55	221.33	1.04	3.58% faster
performance.compilation	383.81	377.9	1.02	1.54% faster
Mean result: 1.002456378048066

Jenkins Console Log
Jenkins Build Stages
Commit hash: 5e38df039f998a983647e228a5501f2cfbeeeac8

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           43 bits physical, 48 bits virtual
CPU(s):                                  256
On-line CPU(s) list:                     0-255
Thread(s) per core:                      2
Core(s) per socket:                      64
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              23
Model:                                   49
Model name:                              AMD EPYC 7742 64-Core Processor
Stepping:                                0
Frequency boost:                         enabled
CPU MHz:                                 1497.185
CPU max MHz:                             3416.0681
CPU min MHz:                             1500.0000
BogoMIPS:                                4491.56
Virtualization:                          AMD-V
L1d cache:                               4 MiB
L1i cache:                               4 MiB
L2 cache:                                64 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-63,128-191
NUMA node1 CPU(s):                       64-127,192-255
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Mitigation; untrained return thunk; SMT enabled with STIBP protection
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Retpolines; IBPB conditional; STIBP always-on; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif v_spec_ctrl umip rdpid overflow_recov succor smca sev sev_es

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[WardBrian]

Even though the functionality is very similar, we prefer to have these in separate files generally