Skip to content

Commit

Permalink
speed up and comment on lyapunov
Browse files Browse the repository at this point in the history
  • Loading branch information
@thorek1
thorek1 committed Jul 31, 2024
1 parent d7db67f commit 87c47a5
Showing 1 changed file with 100 additions and 49 deletions.
149 changes: 100 additions & 49 deletions src/lyapunov.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,16 @@
# Available algorithms:
# :doubling - fast and precise
# :lyapunov - fast for small matrices and precise
# :bicgstab - less precise
# :gmres - less precise
# :iterative - slow and precise
# :speedmapping - slow and very precise

function solve_lyapunov_equation(A::AbstractMatrix{Float64},
C::AbstractMatrix{Float64};
lyapunov_algorithm::Symbol = :doubling,
tol::AbstractFloat = 1e-12)
tol::AbstractFloat = 1e-12,
verbose::Bool = false)
if A isa AbstractSparseMatrix
if length(A.nzval) / length(A) > .1 || lyapunov_algorithm == :lyapunov
A = collect(A)
Expand All @@ -18,7 +27,13 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64},
end
end

solve_lyapunov_equation(A, C, Val(lyapunov_algorithm), tol = tol)
X, solved, i, reached_tol = solve_lyapunov_equation(A, C, Val(lyapunov_algorithm), tol = tol)

if verbose
println("Lyapunov equation - converged to tol $tol: $solved; iterations: $i; reached tol: $reached_tol; algorithm: $lyapunov_algorithm")
end

return X, solved
end

function rrule(::typeof(solve_lyapunov_equation),
Expand Down Expand Up @@ -79,40 +94,45 @@ end


function solve_lyapunov_equation(A::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
C::DenseMatrix{Float64},
C::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
::Val{:lyapunov};
tol::AbstractFloat = 1e-12)
𝐂 = MatrixEquations.lyapd(A, C)

solved = isapprox(𝐂, A * 𝐂 * A' + C, rtol = tol)
𝐂¹ = A * 𝐂 * A' + C

reached_tol = β„’.norm(𝐂¹ - 𝐂) / max(β„’.norm(𝐂), β„’.norm(𝐂¹))

return 𝐂, solved # return info on convergence
return 𝐂, reached_tol < tol, 0, reached_tol # return info on convergence
end



function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64},
C::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
::Val{:doubling};
tol::Float64 = 1e-14)
tol::Float64 = 1e-12)
𝐂 = copy(C)
𝐀 = copy(A)
CA = collect(𝐀)
𝐂A = collect(𝐀)
𝐂¹ = copy(C)

max_iter = 500


iters = max_iter

for i in 1:max_iter
# 𝐂¹ .= 𝐀 * 𝐂 * 𝐀' + 𝐂
mul!(CA, 𝐂, 𝐀')
mul!(𝐂¹, 𝐀, CA, 1, 1)
mul!(𝐂A, 𝐂, 𝐀')
mul!(𝐂¹, 𝐀, 𝐂A, 1, 1)

𝐀 *= 𝐀

droptol!(𝐀, eps())

if i > 10 && i % 2 == 0
if isapprox(𝐂¹, 𝐂, rtol = tol)
iters = i
break
end
end
Expand All @@ -121,55 +141,74 @@ function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64},
# 𝐂 = 𝐂¹
end

solved = isapprox(𝐂, A * 𝐂 * A' + C, rtol = tol)
β„’.mul!(𝐂A, 𝐂, A')
β„’.mul!(𝐂¹, A, 𝐂A)
β„’.axpy!(1, C, 𝐂¹)

denom = max(β„’.norm(𝐂), β„’.norm(𝐂¹))

β„’.axpy!(-1, 𝐂, 𝐂¹)

return 𝐂, solved # return info on convergence
reached_tol = β„’.norm(𝐂¹) / denom

return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence
end




function solve_lyapunov_equation( A::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
C::DenseMatrix{Float64},
C::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
::Val{:doubling};
tol::Float64 = 1e-14)
tol::Float64 = 1e-12)
𝐂 = copy(C)
𝐂¹ = copy(C)
𝐀 = copy(A)

CA = similar(𝐀)
𝐂A = similar(𝐀)
𝐀² = similar(𝐀)

max_iter = 500


iters = max_iter

for i in 1:max_iter
mul!(CA, 𝐂, 𝐀')
mul!(𝐂¹, 𝐀, CA, 1, 1)
mul!(𝐂A, 𝐂, 𝐀')
mul!(𝐂¹, 𝐀, 𝐂A, 1, 1)

mul!(𝐀², 𝐀, 𝐀)
copyto!(𝐀, 𝐀²)

if i > 10 && i % 2 == 0
if isapprox(𝐂¹, 𝐂, rtol = tol)
iters = i
break
end
end

copyto!(𝐂, 𝐂¹)
end

solved = isapprox(𝐂, 𝐀 * 𝐂 * 𝐀' + C, rtol = tol)
β„’.mul!(𝐂A, 𝐂, A')
β„’.mul!(𝐂¹, A, 𝐂A)
β„’.axpy!(1, C, 𝐂¹)

denom = max(β„’.norm(𝐂), β„’.norm(𝐂¹))

β„’.axpy!(-1, 𝐂, 𝐂¹)

return 𝐂, solved # return info on convergence
reached_tol = β„’.norm(𝐂¹) / denom

return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence
end




function solve_lyapunov_equation(A::AbstractMatrix{Float64},
C::DenseMatrix{Float64},
C::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
::Val{:bicgstab};
tol::Float64 = 1e-12)
tol::Float64 = 1e-8)

tmpΜ„ = similar(C)
𝐗 = similar(C)
Expand All @@ -187,16 +226,23 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64},

copyto!(𝐗, 𝐂)

solved = info.solved
β„’.mul!(tmpΜ„, A, 𝐗 * A')
β„’.axpy!(1, C, tmpΜ„)

denom = max(β„’.norm(𝐗), β„’.norm(tmpΜ„))

β„’.axpy!(-1, 𝐗, tmpΜ„)

return 𝐗, solved, solved # return info on convergence
reached_tol = β„’.norm(tmpΜ„) / denom

return 𝐗, reached_tol < tol, info.niter, reached_tol
end


function solve_lyapunov_equation(A::AbstractMatrix{Float64},
C::DenseMatrix{Float64},
C::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
::Val{:gmres};
tol::Float64 = 1e-12)
tol::Float64 = 1e-8)

tmpΜ„ = similar(C)
𝐗 = similar(C)
Expand All @@ -216,30 +262,40 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64},

copyto!(𝐗, 𝐂)

solved = info.solved
β„’.mul!(tmpΜ„, A, 𝐗 * A')
β„’.axpy!(1, C, tmpΜ„)

denom = max(β„’.norm(𝐗), β„’.norm(tmpΜ„))

β„’.axpy!(-1, 𝐗, tmpΜ„)

return 𝐗, solved # return info on convergence
reached_tol = β„’.norm(tmpΜ„) / denom

return 𝐗, reached_tol < tol, info.niter, reached_tol
end


function solve_lyapunov_equation(A::AbstractMatrix{Float64},
C::DenseMatrix{Float64},
C::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
::Val{:iterative};
tol::AbstractFloat = 1e-14)
tol::AbstractFloat = 1e-12)

𝐂 = copy(C)
𝐂¹ = copy(C)
𝐂A = copy(C)

max_iter = 10000

iters = max_iter

for i in 1:max_iter
β„’.mul!(𝐂A, 𝐂, A')
β„’.mul!(𝐂¹, A, 𝐂A)
β„’.axpy!(1, C, 𝐂¹)

if i % 10 == 0
if isapprox(𝐂¹, 𝐂, rtol = tol)
iters = i
break
end
end
Expand All @@ -251,33 +307,28 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64},
β„’.mul!(𝐂¹, A, 𝐂A)
β„’.axpy!(1, C, 𝐂¹)

solved = isapprox(𝐂¹, 𝐂, rtol = tol)
denom = max(β„’.norm(𝐂), β„’.norm(𝐂¹))

β„’.axpy!(-1, 𝐂, 𝐂¹)

return 𝐂, solved # return info on convergence
reached_tol = β„’.norm(𝐂¹) / denom

return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence
end


function solve_lyapunov_equation(A::AbstractMatrix{Float64},
C::DenseMatrix{Float64},
::Val{:speedmapping};
tol::AbstractFloat = 1e-12)

if !(C isa DenseMatrix)
C = collect(C)
end

CA = similar(C)
C::Union{β„’.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
::Val{:speedmapping};
tol::AbstractFloat = 1e-12)
𝐂A = similar(C)

soll = speedmapping(C;
m! = (X, x) -> begin
β„’.mul!(CA, x, A')
β„’.mul!(X, A, CA)
β„’.mul!(𝐂A, x, A')
β„’.mul!(X, A, 𝐂A)
β„’.axpy!(1, C, X)
end, stabilize = false, maps_limit = 10000, tol = tol)

𝐂 = soll.minimizer

solved = soll.converged
end, stabilize = false, maps_limit = 1000, tol = tol)

return 𝐂, solved
return soll.minimizer, soll.converged, soll.maps, soll.norm_βˆ‡
end

0 comments on commit 87c47a5

Please sign in to comment.