Giter Club home page Giter Club logo

periodicschurdecompositions.jl's People

Contributors

ralphas avatar

Stargazers

avatar avatar avatar avatar avatar avatar avatar avatar

Watchers

avatar avatar

periodicschurdecompositions.jl's Issues

ordschur! fails on BigFloat data

I tried to use pschur and ordschur! on extended precision data. The following code fails on BigFloat data:

using PeriodicSchurDecompositions
using JLD
hpd = load("test.jld","hpd");
hpds = [BigFloat.(hpd[i]) for i in 1:200];
PSFs = PeriodicSchurDecompositions.pschur(hpds,:L);
select = abs.(PSFs.values) .< 1;
PeriodicSchurDecompositions.ordschur!(PSFs,select);

ERROR: MethodError: no method matching rmul!(::SubArray{BigFloat, 2, Matrix{BigFloat}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}, ::LinearAlgebra.HessenbergQ{BigFloat, Matrix{BigFloat}, Vector{BigFloat}, false})
Closest candidates are:
  rmul!(::StridedMatrix{T} where T, ::LinearAlgebra.UpperTriangular{<:Any, <:LinearAlgebra.Adjoint}) at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\triangular.jl:1108
  rmul!(::StridedMatrix{T} where T, ::LinearAlgebra.UnitUpperTriangular{<:Any, <:LinearAlgebra.Adjoint}) at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\triangular.jl:1126
  rmul!(::StridedMatrix{T} where T, ::LinearAlgebra.LowerTriangular{<:Any, <:LinearAlgebra.Adjoint}) at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\triangular.jl:1144
  ...
Stacktrace:
 [1] _swapadjqr!(T1::Matrix{BigFloat}, Ts::Vector{Matrix{BigFloat}}, Zs::Vector{Matrix{BigFloat}}, i1::Int64, p1::Int64, p2::Int64; sylcheck::Bool)
   @ PeriodicSchurDecompositions C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\sylswap.jl:138
 [2] _swapadjqr!
   @ C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\sylswap.jl:11 [inlined]
 [3] _swapschur!
   @ C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\rordschur.jl:253 [inlined]
 [4] _moveblock!(P::PeriodicSchur{BigFloat, Matrix{BigFloat}, Matrix{BigFloat}, Matrix{BigFloat}}, jsrc::Int64, jdest::Int64, wantZ::Bool, Q::Vector{Matrix{BigFloat}})
   @ PeriodicSchurDecompositions C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\rordschur.jl:181
 [5] ordschur!(P::PeriodicSchur{BigFloat, Matrix{BigFloat}, Matrix{BigFloat}, Matrix{BigFloat}}, select::BitVector; wantZ::Bool, Z::Nothing)
   @ PeriodicSchurDecompositions C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\rordschur.jl:98
 [6] ordschur!(P::PeriodicSchur{BigFloat, Matrix{BigFloat}, Matrix{BigFloat}, Matrix{BigFloat}}, select::BitVector)
   @ PeriodicSchurDecompositions C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\rordschur.jl:6
 [7] top-level scope
   @ REPL[26]:1

Loading PSD in Julia 1.9 produces strange messages

I tried recently to install all dependencies of the PeriodicSystems with Julia 1.9.0. When loading PSD I got the following messages:

julia> using PeriodicSchurDecompositions
 │ Package PeriodicSchurDecompositions not found, but a package named PeriodicSchurDecompositions is available from a
 │ registry.
 │ Install package?
 │   (@v1.9) pkg> add PeriodicSchurDecompositions
 └ (y/n/o) [y]:
   Resolving package versions...
    Updating `C:\Users\Andreas\.julia\environments\v1.9\Project.toml`
  [e5aedecb] + PeriodicSchurDecompositions v0.1.3
    Updating `C:\Users\Andreas\.julia\environments\v1.9\Manifest.toml`
⌅ [4c555306] ↓ ArrayLayouts v1.0.4 ⇒ v0.8.18
⌃ [aae01518] ↓ BandedMatrices v0.17.24 ⇒ v0.17.18
⌃ [8e7c35d0] ↓ BlockArrays v0.16.28 ⇒ v0.16.25
⌅ [ffab5731] ↓ BlockBandedMatrices v0.12.0 ⇒ v0.11.10
⌅ [1a297f60] ↓ FillArrays v1.1.0 ⇒ v0.13.11
⌃ [4858937d] ↓ InfiniteArrays v0.12.13 ⇒ v0.12.12
⌃ [cde9dba0] ↓ InfiniteLinearAlgebra v0.6.18 ⇒ v0.6.16
⌅ [5078a376] ↓ LazyArrays v1.1.1 ⇒ v0.22.18
⌅ [a3b82374] ↓ MatrixFactorizations v1.0.0 ⇒ v0.9.6
  [e5aedecb] + PeriodicSchurDecompositions v0.1.3
        Info Packages marked with ⌃ and ⌅ have new versions available, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`
Precompiling project...
  ✓ FillArrays
  ✓ ArrayLayouts
  ✓ LowRankApprox
  ✓ FastTransforms
  ✓ Distributions
  ✓ Optim
  ✓ MatrixFactorizations
  ✓ BlockArrays
  ✓ BandedMatrices
  ✓ Distributions → DistributionsChainRulesCoreExt
  ✓ Symbolics
  ✓ DiffEqBase → DiffEqBaseDistributionsExt
  ✓ LazyArrays
  ✓ InfiniteArrays
  ✓ ArrayInterface → ArrayInterfaceBandedMatricesExt
  ✓ BlockBandedMatrices
  ✓ ArrayInterface → ArrayInterfaceBlockBandedMatricesExt
  ✓ LazyBandedMatrices
  ✓ SemiseparableMatrices
  ✓ InfiniteLinearAlgebra
  ✓ ApproxFunBase
  ✓ ApproxFunFourier
  ✓ ApproxFunOrthogonalPolynomials
  ✓ ApproxFunSingularities
  ✓ ApproxFun
  26 dependencies successfully precompiled in 79 seconds. 198 already precompiled.
  25 dependencies precompiled but different versions are currently loaded. Restart julia to access the new versions
┌ Warning: C:\Users\Andreas\.julia\packages\FillArrays\KpKaA\src\oneelement.jl no longer exists, deleted all methods
└ @ Revise C:\Users\Andreas\.julia\packages\Revise\Kw99o\src\packagedef.jl:663
┌ Error: Failed to revise C:\Users\Andreas\.julia\packages\BandedMatrices\qzqxQ\src\symbanded\bandedeigen.jl
│   exception =
│    invalid redefinition of constant BandedEigenvectors
│    Stacktrace:
│     [1] top-level scope
│       @ C:\Users\Andreas\.julia\packages\BandedMatrices\qzqxQ\src\symbanded\bandedeigen.jl:2
│    Revise evaluation error at C:\Users\Andreas\.julia\packages\BandedMatrices\qzqxQ\src\symbanded\bandedeigen.jl:2
│
└ @ Revise C:\Users\Andreas\.julia\packages\Revise\Kw99o\src\packagedef.jl:719
┌ Warning: The running code does not match the saved version for the following files:
│
│   C:\Users\Andreas\.julia\packages\BandedMatrices\qzqxQ\src\symbanded\bandedeigen.jl
│
│ If the error was due to evaluation order, it can sometimes be resolved by calling `Revise.retry()`.
│ Use Revise.errors() to report errors again. Only the first error in each file is shown.
│ Your prompt color may be yellow until the errors are resolved.
└ @ Revise C:\Users\Andreas\.julia\packages\Revise\Kw99o\src\packagedef.jl:829
[ Info: Precompiling PeriodicSchurDecompositions [e5aedecb-f6c0-4c91-b6ff-fbae4296f459]

I wonder if I have to undertake any action. After a new start of Julia using PeriodicSchurDecompositions works, but I wonder which version of the above packages are used.

Registration

I am using the eigenvalue reordering features for a real product to solve continuous-time periodic Riccati equations, as an alternative method to "fast" structure exploiting reductions used in conjunction with multiple-shooting techniques . This functionality is planned to be included in the release 0.4.2 of PeriodicSystems. So, I wonder if the registration of this package is planned in the near future, to allow an easier inclusion of this package.

A second issue is related to supporting the generalized periodic Schur decomposition in the real case. This decomposition would be very useful in both continuous- and discrete-time contexts. Is it is planned to implement the real version of the generalized decomposition?

Here are some technical details. A more efficient solution of the continous-time periodic Riccati equation, would be possible by exploiting features which are intrinsic to the periodic QZ algorithm. The solution technique involves the operation on a left-oriented product of real state transition matrices, formed for a linear periodic system with time-dependent periodic Hamiltonian coefficient matrix. The product S itself is a symplectic matrix, with a dichotomic eigenvalue structure (half of eigenvalues are in the unit circle, half of them outside). In my current implementation, I am moving the "stable" eigenvalues (i.e., those in the unit circle) in the leading positions and operate with the corresponding stable Schur vectors. A known fact in this context is that the "unstable" eigenvalues have the natural tendency to result in the leading positions, so the reordering phase may involve a substantial reordering effort. The methods for constant discrete-time systems try to exploit this feature by working with a symplectic matrix pair (A,E), for which the natural appearance can be exploited by simply interchanging the matrices, i.e., working with (E,A). For the periodic case, this would involve to work with inv(S) (i.e., with the product of inverses), which is, in principle, possible using the generalized decomposition. I am fully exploiting this aspect in the implemetation of "fast" method.
The generalized periodic Schur decomposition in the real case would be also very useful to solve discrete-time periodic Riccati equations.

PS. I guess, I could manage using the available complex decomposition, but this would not be a very elegant and efficient solution.

Erroneous eigenvalue computation

The following eigenvalue computation is part of solving periodic Riccati equations. The undelying periodic matrix hpd, is a 200 vector of 8x8 state transition matrices, whose left product is a simplectic matrix. Thus, the eigenvalues are symmetrical with respect to the unit circle.

The first computation has been performed with the SLICOT wrapper pschur available in the PeriodicSystems and produces the following eigenvalues:

using PeriodicSystems
using JLD
hpd = load("test.jld","hpd")
S, Z, ev, = PeriodicSystems.pschur(hpd); ev



8-element Vector{ComplexF64}:
   9.473709584260244e71 + 6.144171599316869e71im
   9.473709584260244e71 - 6.144171599316869e71im
   7.833585365370729e44 + 0.0im
    4.644872743033464e6 + 0.0im
  2.1529115119486817e-7 + 0.0im
 1.2765546724226738e-45 + 0.0im
  7.430242860588899e-73 + 4.8188818491875555e-73im
  7.430242860588899e-73 - 4.8188818491875555e-73im

It can be easily checked the symmetry of eigenvalues

julia> [1. ./ev ev]
8×2 Matrix{ComplexF64}:
 7.43024e-73-4.81888e-73im   9.47371e71+6.14417e71im
 7.43024e-73+4.81888e-73im   9.47371e71-6.14417e71im
 1.27655e-45-0.0im           7.83359e44+0.0im
  2.15291e-7-0.0im            4.64487e6+0.0im
   4.64487e6-0.0im           2.15291e-7+0.0im
  7.83359e44-0.0im          1.27655e-45+0.0im
  9.47371e71-6.14417e71im   7.43024e-73+4.81888e-73im
  9.47371e71+6.14417e71im   7.43024e-73-4.81888e-73im

The same computation performed with the PeriodicSchurDecompositions produces wrong eigenvalues:

using PeriodicSchurDecompositions
using JLD
hpd = load("test.jld","hpd")
PSF = PeriodicSchurDecompositions.pschur(hpd,:L); 

julia> PSF.values
8-element Vector{ComplexF64}:
   9.473709584260287e71 + 6.144171599318682e71im
   9.473709584260287e71 - 6.144171599318682e71im
                    0.0 + 0.0im
   3.184085626273608e45 + 0.0im
  -9.313225746154785e-9 + 0.0im
   4.0418154850413655e6 + 0.0im
 2.3165926597911904e-45 + 0.0im
                    0.0 + 0.0im

The next computation would be the reordering of eigenvalues, but this is probably another story.
I am attaching the data in a separate file.
test.zip

Eigenvectors

Hi,

Is there a simple / generic way to retrieve the eigenvectors?

Extended decompositions (time-varying matrix size)

As described in #2 handling time-varying matrix size would be a valuable addition to this package. On looking into it the required algorithms seem straight-forward, so it may get done in the coming months. I'm adding an issue to track (eventual) progress.

failures with unit period

Hi,

I think there is a corner case that is missing

partial_pschur([rand(100,100)+I],10) #errors
partial_pschur([rand(100,100)+I,rand(100,100)+I],10) #does not

Accuracy loss for small magnitude eigenvalues

This is basically a similar example to that in #4, but for a symplectic product of K = 300 8x8 matrices (instead K = 200) (the underlying Riccati equation is the same, but the number of discretization steps K is different, i.e. 300 vs. 200).

The first computation, to serve as reference, has been performed with the SLICOT wrapper pschur available in the PeriodicSystems and produces the following eigenvalues and their reciprocals:

using PeriodicSystems
using JLD
hpd = load("test1.jld","hpd");
S, Z, ev, sind, = PeriodicSystems.pschur(hpd); 
[sort(ev,by=abs) sort(1 ./ev,by=abs)] 

julia> [sort(ev,by=abs) sort(1 ./ev,by=abs)]
8×2 Matrix{ComplexF64}:
 -1.92309e-56+4.89783e-56im  -1.92309e-56-4.89783e-56im
 -1.92309e-56-4.89783e-56im  -1.92309e-56+4.89783e-56im
  2.43966e-47+0.0im           2.43966e-47-0.0im
   2.21356e-7+0.0im            2.21356e-7-0.0im
    4.51761e6+0.0im             4.51761e6-0.0im
   4.09892e46+0.0im            4.09892e46-0.0im
  -6.94581e54+1.769e55im      -6.94581e54-1.769e55im
  -6.94581e54-1.769e55im      -6.94581e54+1.769e55im

which clearly illustrates the symplectic nature of the eigenproblem.

The same computation performed with the PeriodicSchurDecompositions produces slightly inaccurate eigenvalues, which do not fullfil the symplectic requirement, as can be seen below:

using PeriodicSchurDecompositions
using JLD
hpd = load("test1.jld","hpd");
PSF = PeriodicSchurDecompositions.pschur(hpd,:L);
ev1 = PSF.values;
[sort(ev1,by=abs) sort(1 ./ev1,by=abs)]


julia> [sort(ev1,by=abs) sort(1 ./ev1,by=abs)]
8×2 Matrix{ComplexF64}:
  1.07037e-56+0.0im       -1.92309e-56-4.89783e-56im
 -1.65442e-56+0.0im       -1.92309e-56+4.89783e-56im
  2.36581e-48-0.0im        2.43966e-47-0.0im
   2.21356e-7-0.0im         2.21356e-7-0.0im
    4.51761e6+0.0im          4.51761e6+0.0im
   4.09892e46+0.0im         4.22688e47+0.0im
  -6.94581e54+1.769e55im   -6.04442e55-0.0im
  -6.94581e54-1.769e55im    9.34258e55-0.0im

This loss of accuracy in the computed small size eigenvalues prevents the success of some tests for solving periodic Riccati differential equations. The computations for K = 100 are OK on Windows and Linux with Julia 1.7, but fail for K = 200 and K = 300 on both Windows and Linux.

I am attaching the data in the zip-file
test1.zip

Vector of matrices vs. 3-dimensional array

What is the reason preferring the representation of periodic matrices via vectors of matrices over 3-dimensional arrays (as in the original SLICOT implementation)?

The vector of matrices based implementation would in principle allow working with time-varying dimensions. Do you intend to support at certain later time time-varying dimensions?

matrixfree partial_gpschur

Hi,

Thanks a lot for your work!

Would it make sense to have a Matrix-free generalized periodic Schur decomposition? How would one provide Bₚ⁻¹? Just Bₚ?

Or maybe it is better to use partial_pschur where the user implement x -> B \ A*x?

Thanks for your input

More generic interface for pschur!

The interface of pschur! is currently defined as

    pschur!(A::Vector{<:StridedMatrix}, S::Vector{Bool}, lr::Symbol) -> F::GeneralizedPeriodicSchur

Computes a generalized periodic Schur decomposition of a series of general square matrices
with left (`lr=:L`) or right (`lr=:R`) orientation.
Entries in `S` must be `false` for the corresponding entries in `A` which are effectively
inverted, and `true` for the rest. Currently `Sⱼ` must be `true` for the leftmost term.

I wonder if it would be possible to get rid of the restriction that

Currently Sⱼmust betrue for the leftmost term.

I can imagine that removing this restriction is not trivial and can even lead to the redefinition of the GeneralizedPeriodicSchur object
and even to new computational approaches. However, for user's convenience this would be (in my opinion) very useful.

I implemented myself a similar interface routine to the SLICOT wrappers, without such a restriction (this was relatively easy, because the SLICOT wrapper to MB03BD supports exactly this). This allowed me, for example, to easily switch in the recently implemented periodic Riccati solvers from a quotient-product form to its inverse (where the Schur element is placed into a matrix which enters as an inverse factor).

I am watching with much interest your developments, because I am seriously concerned of the difficulties to generate correct binary libraries on a range of compilers and OSs for the SLICOT library. Although the SLICOT tools are very reliable and fast, I am pondering on the long range to get rid of SLICOT binding, provided alternative tools will be available in your package. I already implemented an alternative SLICOT-free interface for the continuous-time Riccati solvers, using the interface

pschur!(A::Vector{S<:StridedMatrix}, lr::Symbol=:R) -> F::PeriodicSchur

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.