From e11e62e7f730d88afe73f591394650734c1009f3 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Wed, 18 Dec 2024 20:05:48 +0100 Subject: [PATCH] docs for lyap and sylv. no more iterative. rename bartels_stewart --- src/MacroModelling.jl | 2 +- src/algorithms/lyapunov.jl | 81 +++++++++++++------------- src/algorithms/sylvester.jl | 111 ++++++++++++++++++------------------ src/common_docstrings.jl | 2 + src/get_functions.jl | 15 ++++- src/perturbation.jl | 2 +- 6 files changed, 113 insertions(+), 100 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 7e27e931..60ea2727 100644 --- 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 index d9ae8ba0..09ab1225 100644 --- 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 index 5576e6b0..a385610b 100644 --- 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,7 +67,7 @@ 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) @@ -74,7 +75,7 @@ function solve_sylvester_equation(A::M, 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 index fb7fc231..6ecf5dce 100644 --- 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`" \ No newline at end of file diff --git a/src/get_functions.jl b/src/get_functions.jl index ca1acbe2..1176028c 100644 --- 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 index 0e2686eb..7ce85fdd 100644 --- 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()