docs for lyap and sylv. no more iterative. rename bartels_stewart

thorek1 · Dec 18, 2024 · e11e62e · e11e62e
1 parent 3e9fcb1
commit e11e62e
Show file tree

Hide file tree

Showing 6 changed files with 113 additions and 100 deletions.
diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl
@@ -7406,7 +7406,7 @@ function get_relevant_steady_state_and_state_update(::Val{:first_order},
                                                     𝓂::ℳ, 
                                                     tol::AbstractFloat; 
                                                     quadratic_matrix_equation_solver::Symbol = :schur, 
-                                                    sylvester_algorithm::Symbol = :gmres, 
+                                                    sylvester_algorithm::Symbol = :bicgstab, 
                                                     verbose::Bool = false,
                                                     timer::TimerOutput = TimerOutput())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{Float64}}, Bool} where S <: Real
     SS_and_pars, (solution_error, iters) = get_NSSS_and_parameters(𝓂, parameter_values, tol = tol, timer = timer, verbose = verbose)

diff --git a/src/algorithms/lyapunov.jl b/src/algorithms/lyapunov.jl
@@ -1,8 +1,9 @@
 # Available algorithms: 
 # :doubling     - fast and precise
-# :lyapunov     - fast for small matrices and precise, dense matrices only
+# :bartels_stewart     - fast for small matrices and precise, dense matrices only
 # :bicgstab     - less precise
 # :gmres        - less precise
+
 # :iterative    - slow and precise
 # :speedmapping - slow and very precise
 
@@ -17,7 +18,7 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64},
     @timeit_debug timer "Solve lyapunov equation" begin
     @timeit_debug timer "Choose matrix formats" begin
 
-    if lyapunov_algorithm ≠ :lyapunov
+    if lyapunov_algorithm ≠ :bartels_stewart
         A = choose_matrix_format(A)
     end
 
@@ -54,7 +55,7 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64},
             C = collect(C)
 
             X, i, reached_tol = solve_lyapunov_equation(A, C, 
-                                                                Val(:lyapunov), 
+                                                                Val(:bartels_stewart), 
                                                                 # tol = tol, 
                                                                 timer = timer)
 
@@ -145,7 +146,7 @@ end
 
 function solve_lyapunov_equation(A::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
                                 C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
-                                ::Val{:lyapunov};
+                                ::Val{:bartels_stewart};
                                 tol::AbstractFloat = 1e-12,
                                 timer::TimerOutput = TimerOutput())
     𝐂 = try 
@@ -469,54 +470,54 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64},
 end
 
 
-function solve_lyapunov_equation(A::AbstractMatrix{Float64},
-                                C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
-                                ::Val{:iterative};
-                                tol::AbstractFloat = 1e-14,
-                                timer::TimerOutput = TimerOutput())
-    𝐂  = copy(C)
-    𝐂¹ = copy(C)
-    𝐂A = copy(C)
+# function solve_lyapunov_equation(A::AbstractMatrix{Float64},
+#                                 C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
+#                                 ::Val{:iterative};
+#                                 tol::AbstractFloat = 1e-14,
+#                                 timer::TimerOutput = TimerOutput())
+#     𝐂  = copy(C)
+#     𝐂¹ = copy(C)
+#     𝐂A = copy(C)
 
-    max_iter = 10000
+#     max_iter = 10000
 
-    iters = max_iter
+#     iters = max_iter
 
-    for i in 1:max_iter
-        ℒ.mul!(𝐂A, 𝐂, A')
-        ℒ.mul!(𝐂¹, A, 𝐂A)
-        ℒ.axpy!(1, C, 𝐂¹)
+#     for i in 1:max_iter
+#         ℒ.mul!(𝐂A, 𝐂, A')
+#         ℒ.mul!(𝐂¹, A, 𝐂A)
+#         ℒ.axpy!(1, C, 𝐂¹)
 
-        if i % 10 == 0
-            normdiff = ℒ.norm(𝐂¹ - 𝐂)
-            if !isfinite(normdiff) || normdiff / max(ℒ.norm(𝐂), ℒ.norm(𝐂¹)) < tol
-            # if isapprox(𝐂¹, 𝐂, rtol = tol)
-                iters = i
-                break
-            end
-        end
+#         if i % 10 == 0
+#             normdiff = ℒ.norm(𝐂¹ - 𝐂)
+#             if !isfinite(normdiff) || normdiff / max(ℒ.norm(𝐂), ℒ.norm(𝐂¹)) < tol
+#             # if isapprox(𝐂¹, 𝐂, rtol = tol)
+#                 iters = i
+#                 break
+#             end
+#         end
 
-        copyto!(𝐂, 𝐂¹)
-    end
+#         copyto!(𝐂, 𝐂¹)
+#     end
 
-    # ℒ.mul!(𝐂A, 𝐂, A')
-    # ℒ.mul!(𝐂¹, A, 𝐂A)
-    # ℒ.axpy!(1, C, 𝐂¹)
+#     # ℒ.mul!(𝐂A, 𝐂, A')
+#     # ℒ.mul!(𝐂¹, A, 𝐂A)
+#     # ℒ.axpy!(1, C, 𝐂¹)
 
-    # denom = max(ℒ.norm(𝐂), ℒ.norm(𝐂¹))
+#     # denom = max(ℒ.norm(𝐂), ℒ.norm(𝐂¹))
 
-    # ℒ.axpy!(-1, 𝐂, 𝐂¹)
+#     # ℒ.axpy!(-1, 𝐂, 𝐂¹)
 
-    # reached_tol = denom == 0 ? 0.0 : ℒ.norm(𝐂¹) / denom
+#     # reached_tol = denom == 0 ? 0.0 : ℒ.norm(𝐂¹) / denom
 
-    reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂)
+#     reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂)
 
-    # if reached_tol > tol
-    #     println("Lyapunov: iterative $reached_tol")
-    # end
+#     # if reached_tol > tol
+#     #     println("Lyapunov: iterative $reached_tol")
+#     # end
 
-    return 𝐂, iters, reached_tol # return info on convergence
-end
+#     return 𝐂, iters, reached_tol # return info on convergence
+# end
 
 
 # function solve_lyapunov_equation(A::AbstractMatrix{Float64},

diff --git a/src/algorithms/sylvester.jl b/src/algorithms/sylvester.jl
@@ -1,9 +1,10 @@
 # Available algorithms: 
 # :doubling     - fast, expensive part: B^2
-# :sylvester    - fast, dense matrices only
+# :bartels_stewart    - fast, dense matrices only
 # :bicgstab     - fastest for large problems, might not reach desired precision, warm start not always helpful
 # :dqgmres      - fastest for large problems, might not reach desired precision, stable path with warm start
 # :gmres      - fastest for large problems, might not reach desired precision, can be effective, not efficient
+
 # :iterative    - slow
 # :speedmapping - slow
 
@@ -20,13 +21,13 @@ function solve_sylvester_equation(A::M,
                                     verbose::Bool = false) where {M <: AbstractMatrix{Float64}, N <: AbstractMatrix{Float64}, O <: AbstractMatrix{Float64}}
     @timeit_debug timer "Choose matrix formats" begin
 
-    if sylvester_algorithm == :sylvester
+    if sylvester_algorithm == :bartels_stewart
         b = collect(B)
     else
         b = choose_matrix_format(B)# |> collect
     end
 
-    if sylvester_algorithm ∈ [:bicgstab, :gmres, :sylvester]
+    if sylvester_algorithm ∈ [:bicgstab, :gmres, :bartels_stewart]
         a = collect(A)
 
         c = collect(C)
@@ -66,15 +67,15 @@ function solve_sylvester_equation(A::M,
         println("Sylvester equation - converged to tol $tol: $(reached_tol < tol); iterations: $i; reached tol: $reached_tol; algorithm: $sylvester_algorithm")
     end
 
-    if !(reached_tol < tol) && sylvester_algorithm ≠ :sylvester && length(B) < 5e7 # try sylvester if previous one didn't solve it
+    if !(reached_tol < tol) && sylvester_algorithm ≠ :bartels_stewart && length(B) < 5e7 # try sylvester if previous one didn't solve it
         aa = collect(A)
 
         bb = collect(B)
 
         cc = collect(C)
 
         x, i, reached_tol = solve_sylvester_equation(aa, bb, cc, 
-                                                            Val(:sylvester), 
+                                                            Val(:bartels_stewart), 
                                                             initial_guess = zeros(0,0), 
                                                             # tol = tol, 
                                                             verbose = verbose,
@@ -967,7 +968,7 @@ end
 function solve_sylvester_equation(A::DenseMatrix{Float64},
                                     B::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
                                     C::DenseMatrix{Float64},
-                                    ::Val{:sylvester};
+                                    ::Val{:bartels_stewart};
                                     initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0),
                                     timer::TimerOutput = TimerOutput(),
                                     verbose::Bool = false,
@@ -1337,71 +1338,71 @@ function solve_sylvester_equation(A::DenseMatrix{Float64},
 end
 
 
-function solve_sylvester_equation(A::AbstractMatrix{Float64},
-                                    B::AbstractMatrix{Float64},
-                                    C::AbstractMatrix{Float64},
-                                    ::Val{:iterative};
-                                    initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0),
-                                    timer::TimerOutput = TimerOutput(),
-                                    verbose::Bool = false,
-                                    tol::AbstractFloat = 1e-14)
-    # guess_provided = true
+# function solve_sylvester_equation(A::AbstractMatrix{Float64},
+#                                     B::AbstractMatrix{Float64},
+#                                     C::AbstractMatrix{Float64},
+#                                     ::Val{:iterative};
+#                                     initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0),
+#                                     timer::TimerOutput = TimerOutput(),
+#                                     verbose::Bool = false,
+#                                     tol::AbstractFloat = 1e-14)
+#     # guess_provided = true
 
-    if length(initial_guess) == 0
-        # guess_provided = false
-        initial_guess = zero(C)
-    end
+#     if length(initial_guess) == 0
+#         # guess_provided = false
+#         initial_guess = zero(C)
+#     end
 
-    𝐂  = A * initial_guess * B + C - initial_guess
-    𝐂⁰  = copy(𝐂)
-    # 𝐂  = copy(C)
+#     𝐂  = A * initial_guess * B + C - initial_guess
+#     𝐂⁰  = copy(𝐂)
+#     # 𝐂  = copy(C)
 
-    𝐂¹ = similar(C)
-    𝐂B = similar(C)
+#     𝐂¹ = similar(C)
+#     𝐂B = similar(C)
 
-    max_iter = 10000
+#     max_iter = 10000
 
-    iters = max_iter
+#     iters = max_iter
 
-    for i in 1:max_iter
-        @timeit_debug timer "Update" begin
-        ℒ.mul!(𝐂B, 𝐂, B)
-        ℒ.mul!(𝐂¹, A, 𝐂B)
-        ℒ.axpy!(1, 𝐂⁰, 𝐂¹)
+#     for i in 1:max_iter
+#         @timeit_debug timer "Update" begin
+#         ℒ.mul!(𝐂B, 𝐂, B)
+#         ℒ.mul!(𝐂¹, A, 𝐂B)
+#         ℒ.axpy!(1, 𝐂⁰, 𝐂¹)
 
-        if i % 10 == 0
-            normdiff = ℒ.norm(𝐂¹ - 𝐂)
-            if !isfinite(normdiff) || normdiff / max(ℒ.norm(𝐂), ℒ.norm(𝐂¹)) < tol
-            # if isapprox(𝐂¹, 𝐂, rtol = tol)
-                iters = i
-                break
-            end
-        end
+#         if i % 10 == 0
+#             normdiff = ℒ.norm(𝐂¹ - 𝐂)
+#             if !isfinite(normdiff) || normdiff / max(ℒ.norm(𝐂), ℒ.norm(𝐂¹)) < tol
+#             # if isapprox(𝐂¹, 𝐂, rtol = tol)
+#                 iters = i
+#                 break
+#             end
+#         end
 
-        copyto!(𝐂, 𝐂¹)
-        end # timeit_debug
-    end
+#         copyto!(𝐂, 𝐂¹)
+#         end # timeit_debug
+#     end
 
-    # ℒ.mul!(𝐂B, 𝐂, B)
-    # ℒ.mul!(𝐂¹, A, 𝐂B)
-    # ℒ.axpy!(1, C, 𝐂¹)
+#     # ℒ.mul!(𝐂B, 𝐂, B)
+#     # ℒ.mul!(𝐂¹, A, 𝐂B)
+#     # ℒ.axpy!(1, C, 𝐂¹)
 
-    # denom = max(ℒ.norm(𝐂), ℒ.norm(𝐂¹))
+#     # denom = max(ℒ.norm(𝐂), ℒ.norm(𝐂¹))
 
-    # ℒ.axpy!(-1, 𝐂, 𝐂¹)
+#     # ℒ.axpy!(-1, 𝐂, 𝐂¹)
 
-    # reached_tol = denom == 0 ? 0.0 : ℒ.norm(𝐂¹) / denom
+#     # reached_tol = denom == 0 ? 0.0 : ℒ.norm(𝐂¹) / denom
 
-    𝐂 += initial_guess
+#     𝐂 += initial_guess
 
-    reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C))
+#     reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C))
 
-    # if reached_tol > tol
-    #     println("Sylvester: iterative $reached_tol")
-    # end
+#     # if reached_tol > tol
+#     #     println("Sylvester: iterative $reached_tol")
+#     # end
 
-    return 𝐂, iters, reached_tol # return info on convergence
-end
+#     return 𝐂, iters, reached_tol # return info on convergence
+# end
 
 
 # function solve_sylvester_equation(A::AbstractMatrix{Float64},

diff --git a/src/common_docstrings.jl b/src/common_docstrings.jl
@@ -19,3 +19,5 @@ const PARAMETER_DERIVATIVES = "`parameter_derivatives` [Default: :all]: paramete
 const DATA = "`data` [Type: `KeyedArray`]: data matrix with variables (`String` or `Symbol`) in rows and time in columns"
 const SMOOTH = "`smooth` [Default: `true`, Type: `Bool`]: whether to return smoothed (`true`) or filtered (`false`) shocks/variables. Only works for the Kalman filter. The inversion filter only returns filtered shocks/variables."
 const DATA_IN_LEVELS = "`data_in_levels` [Default: `true`, Type: `Bool`]: indicator whether the data is provided in levels. If `true` the input to the data argument will have the non stochastic steady state substracted."
+const LYAPUNOV = "`lyapunov_algorithm` [Default: `:doubling`, Type: `Symbol`]: algorithm to solve Lyapunov equation (`A * X * A' + C = X`). Available algorithms: `:doubling`, `:bartels_stewart`, `:bicgstab`, `:gmres`"
+const SYLVESTER = "`sylvester_algorithm` [Default: `:doubling`, Type: `Symbol`]: algorithm to solve Sylvester equation (`A * X * B + C = X`). Available algorithms: `:doubling`, `:bartels_stewart`, `:bicgstab`, `:dqgmres`, `:gmres`"
diff --git a/src/get_functions.jl b/src/get_functions.jl
@@ -1974,6 +1974,7 @@ Return the variance decomposition of endogenous variables with regards to the sh
 - $MODEL
 # Keyword Arguments
 - $PARAMETERS
+- $LYAPUNOV
 - $VERBOSE
 
 # Examples
@@ -2086,6 +2087,8 @@ Return the correlations of endogenous variables using the first, pruned second,
 # Keyword Arguments
 - $PARAMETERS
 - $ALGORITHM
+- $LYAPUNOV
+- $SYLVESTER
 - $VERBOSE
 
 # Examples
@@ -2179,6 +2182,8 @@ Return the autocorrelations of endogenous variables using the first, pruned seco
 - `autocorrelation_periods` [Default: `1:5`]: periods for which to return the autocorrelation
 - $PARAMETERS
 - $ALGORITHM
+- $LYAPUNOV
+- $SYLVESTER
 - $VERBOSE
 
 # Examples
@@ -2290,6 +2295,8 @@ Return the first and second moments of endogenous variables using the first, pru
 - $DERIVATIVES
 - $PARAMETER_DERIVATIVES
 - $ALGORITHM
+- $LYAPUNOV
+- $SYLVESTER
 - `dependencies_tol` [Default: `1e-12`, Type: `AbstractFloat`]: tolerance for the effect of a variable on the variable of interest when isolating part of the system for calculating covariance related statistics
 - $VERBOSE
 
@@ -2798,6 +2805,8 @@ Function to use when differentiating model moments with repect to parameters.
 - `autocorrelation` [Default: `Symbol[]`, Type: `Vector{Symbol}`]: if values are provided the function returns the autocorrelation of the mentioned variables
 - `autocorrelation_periods` [Default: `1:5`]: periods for which to return the autocorrelation of the mentioned variables
 - $ALGORITHM
+- $LYAPUNOV
+- $SYLVESTER
 - $VERBOSE
 
 # Examples
@@ -3032,9 +3041,9 @@ function get_loglikelihood(𝓂::ℳ,
     timer::TimerOutput = TimerOutput(),
     verbose::Bool = false)::S where S <: Real
 
-    if algorithm ∈ [:third_order,:pruned_third_order]
-        sylvester_algorithm = :bicgstab
-    end
+    # if algorithm ∈ [:third_order,:pruned_third_order]
+    #     sylvester_algorithm = :bicgstab
+    # end
 
     # TODO: throw error for bounds violations, suggesting this might be due to wrong parameter ordering
     @assert length(parameter_values) == length(𝓂.parameters) "The number of parameter values provided does not match the number of parameters in the model. If this function is used in the context of estimation and not all parameters are estimated, you need to combine the estimated parameters with the other model parameters in one `Vector`. Make sure they have the same order they were declared in the `@parameters` block (check by calling `get_parameters`)."

diff --git a/src/perturbation.jl b/src/perturbation.jl
@@ -255,7 +255,7 @@ function rrule(::typeof(calculate_first_order_solution),
 
         tmp1 = M' * ∂𝐒ᵗ * expand[2]
 
-        ss, solved = solve_sylvester_equation(tmp2, 𝐒̂ᵗ', -tmp1, sylvester_algorithm = :sylvester)
+        ss, solved = solve_sylvester_equation(tmp2, 𝐒̂ᵗ', -tmp1, sylvester_algorithm = :doubling)
 
         if !solved
             NoTangent(), NoTangent(), NoTangent()