Adjoint handling

[dkarrasch]

From what I understand, the adj field is typically computed from type information, like isa(A, Adjoint) and things like that. That means that the value of adj can be computed at compile time. Once the objects are generated, this type information is lost and only a Bool value is kept. This makes it unreachable for the compiler and in particular for multiple dispatch.

I'd propose to consider to move the adj information into a type parameter, for instance like:

struct GeneralizedLyapunovMap{T,TA <: AbstractMatrix{T},TE <: AbstractMatrix{T},Adj} <: LyapunovMatrixEquationsMaps{T}
   A::TA
   E::TE
   disc::Bool
   her::Bool
   adj::Bool
   function GeneralizedLyapunovMap(A::AbstractMatrix{T}, E::AbstractMatrix{T}; disc=false, her=false) where {T}
      n = LinearAlgebra.checksquare(A)
      n == LinearAlgebra.checksquare(E) ||
               throw(DimensionMismatch("E must be a square matrix of dimension $n"))
      return new{T,typeof(A),typeof(E),false}(A, E, disc, her)
   end
   function GeneralizedLyapunovMap(A::Adjoint{T}, E::Adjoint{T}; disc=false, her=false) where {T}
      n = LinearAlgebra.checksquare(A)
      n == LinearAlgebra.checksquare(E) ||
               throw(DimensionMismatch("E must be a square matrix of dimension $n"))
      return new{T,typeof(A),typeof(E),true}(A, E, disc, her)
   end
end

You could then dispatch on types ::GeneralizedLyapunovMap{T,TA,TE,true} and ::GeneralizedLyapunovMap{T,TA,TE,false} if needed. In this specific case, I don't even see why you would need that information, at least not in the meoperators part. Anyway, I don't think it would lead to a major performance increase, but it could help clean up the code a little bit, i.e., have less branches.

[dkarrasch]

Thinking it further, some other "values" could be replaced by types:

abstract type DiscreteOrContinuous end
struct Discrete <: DiscreteOrContinuous end
struct Continuous <: DiscreteOrContinuous end

Now you could store the previous disc as a type parameter instead of as the value of a field. That would remove two nested branches in some cases and lead to short, concise methods. If you want, the branching would then be performed by multiple dispatch.

On a different note, why do you have the struct-constructors and then these smallcased, abbreviated functions like lyapop? The latter look like outer constructors, which could (should?) have the same name like the type you're constructing. Just throwing out ideas...

[andreasvarga]

Thanks a lot for these valuable comments. I am about to release version v2.0 of ME announcing also the transition to LM. My last effort was to improve the performance of some of the basic solvers, so I ignored to check the issues. Sorry! (I wonder if it is any possibility to automatically get the notifications of new issues or PRs?)

I think I need some time to evaluate the implications of your proposals and probably I will delay releasing v2.0 (it was planned for today, before my one week vacancy).

Regarding the layer of functions like lyapop, this has historical reasons. I tried to use the same names to define operators as I used previously, when I employed LinearOperators. So, there is no change in the interface to the previous versions, which I think is good for potential users. Moreover, this allows to use the same operator name for Lyapunov operators defined for one matrix (e.g., L = A*X+X*A') and a pair of matrices (L = A*X*E'+E*X*A'). So, I was able to keep the documentation practically unchanged.

The handling of discrete/continuous issue as you propose, would probably have some implications at the level of basic solvers, which now have even different names (e.g.,lyapc/lyapd).

BTW: I was able to define for all my operators also the corresponding inverse operators, callable via the LinearAlgebra.inv. This allowed to implement the functionrcondest for reverse condition number estimation using a single operator as input parameter (previously I used two operators, where one was assumed to be the inverse of the other). I wonder if it would be appropriate to add to LM functions to compute the L1-norm of operators (exactly or approximately as done in LAPACK).

[andreasvarga]

Regarding the Adjoint issue, at a certain moment I even removed the adj parameter from the definition of Lyapunov operators. This would also be possible for the inverse operators, where this information is used only in the case when the matrices are in Schur/generalized Schur form, since this involves the use of alternative lower level (and faster) solvers which use this information. In the previous version relying on LinearOperators I had separate functions to construct operators for Schur form matrices, so I am happy that it was possible to simplify the definitions by automatically recognizing if the matrices are reduced or not.

[dkarrasch]

I wonder if it is any possibility to automatically get the notifications of new issues or PRs?

I think you need to change the "watch" status at the top of the github page of ME.jl.

The handling of discrete/continuous issue as you propose, would probably have some implications at the level of basic solvers, which now have even different names (e.g., lyapc/lyapd).

Not necessarily. You could simply dispatch on that value, let's say

somefunction(::GeneralizedLyapunovMap{T,TA,TE,true}, args...; kwargs...) where {T,TA,TE}
    # code that assumes
    adj = true
end

somefunction(::GeneralizedLyapunovMap{T,TA,TE,false}, args...; kwargs...) where {T,TA,TE}
    # code that assumes
    adj = false
end

This reduces the amount of branching, since branching would be done by multiple dispatch and not at runtime.

[andreasvarga]

Thanks for the hints.

Since in the meantime I made some modifications (I think the operator part was not effected), I wonder how to handle the PRs which you suggested. Is it possible that some conflicts occur?

[dkarrasch]

Yes, there are some conflicts. I'll resolve them, but then it would be good to merge the PR and continue from there. 😅

[andreasvarga]

Thanks.

[dkarrasch]

So, for the LyapunovMap, I don't see any usage of the L.adj value anywhere in the package. Could this be removed? Same with GeneralizedLyapunovMap. The first usage I see is with InverseLyapunovMap.

[andreasvarga]

I mentioned that for some reason, after initially removing adj from the definitions of LyapunovMap and GeneralizedLyapunovMap, I reintroduced it again. I think I intended to ensure a certain symmetry in the way how I overloaded the LinearAlgebra.inv function for LyapunovMap:

LinearAlgebra.inv(L::LyapunovMap) = invlyapop(L.A; disc=L.disc, her=L.her)
LinearAlgebra.inv(L::InverseLyapunovMap) = lyapop(L.A; disc=L.disc, her=L.her)

When forming the inverse map withinvlyapop, the information on adjoint is stored in adj and L.A contains the non-transposed array, that is, either A.parent or A. So, I now realize that the second line is wrong and should be:

LinearAlgebra.inv(L::InverseLyapunovMap) = lyapop(L.adj ? L.A' : L.A; disc=L.disc, her=L.her)

The same applies for GeneralizedLyapunovMap .

I also think adj can be removed from the LyapunovMap and GeneralizedLyapunovMap structures.

[dkarrasch]

I included that in my ongoing PR.

[andreasvarga]

I will not be able to perform any operation, probably until next Sunday (25.07.2021). I checked that similar changes of LinearAlgebra.inv are also necessary for the Sylvester maps.

However, I am now receiving all notifications.

[dkarrasch]

No worries. Have a good vacation.

[andreasvarga]

I am back. I had a look to your changes. For me they are OK. Should I merge them?

[dkarrasch]

Did you try out some benchmarks?

[andreasvarga]

No. Actually I have no dedicated benchmarks for the operator functions. I intended to perform some tests after merging and pulling the new version to my local dev version. Or there is another way to make tests before merging? Sorry for my ignorance!

[dkarrasch]

No worries. If you use VS Code with the github extension, then you can simply "check out a PR". That will download my branch and add a git branch into your dev folder. You can then start a new Julia session and will run my branch, without committing or merging.

[andreasvarga]

I installed the Pull Request feature, which is a nice feature for cooperative developments!

The two broken tests can be actually fixed with the new implementation of operators. The issue is that the real solvers exploit the real Schur form and work only for real right-hand side (no support for complex right-hand side). Therefore, for complex right-hand side, the Schur structure can be simply ignored and the general solvers can be used.

[andreasvarga]

I think I have troubles with the master repository. Erroneously I created a new branch andreasvarga/master which now contains your corrections and my last modifications (I fixed the broken tests). However, the automatic merging into master is apparently not possible and, to be honest, I have no idea what to do to get rid of this branch and integrate it in the master repository. Also locally things are now confused. I guess, my attempt to perform tests with your PR went wrong. I would appreciate very much your help in fixing these issues.

[dkarrasch]

Alright, so my PR is already merged into master, as well as your change to the release notes. If your further changes are not too difficult, you could delete the MatrixEquations folder in your local dev and then start anew: dev MatrixEquations in pkg mode and off you go with creating new branches, make the changes, open a PR, let CI run and merge if all is well.

It could be that you have changed the remote repository to mine while checking out my branch. In VS Code you can simply choose the current branch in the lower left corner, I don't even use/know the command line commands for doing that.

[andreasvarga]

What about deleting everything locally and clone a new repository from the current master?

Daniel Karrasch @.***> schrieb am Mo., 26. Juli 2021, 20:26:

Alright, so my PR is already merged into master, as well as your change to the release notes. If your further changes are not too difficult, you could delete the MatrixEquations folder in your local dev and then start anew: dev MatrixEquations in pkg mode and off you go with creating new branches, make the changes, open a PR, let CI run and merge if all is well.

It could be that you have changed the remote repository to mine while checking out my branch. In VS Code you can simply choose the current branch in the lower left corner, I don't even use/know the command line commands for doing that.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/andreasvarga/MatrixEquations.jl/issues/15#issuecomment-886928252, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEAY5DX5TE3F5F2ZWILTZWSDRANCNFSM5AEFHO6A .

[dkarrasch]

Exactly that's what's happening when you delete locally and dev MatrixEquations again.

[andreasvarga]

Finally, I recovered the desired state and also fixed an initialization issue in gsylvs! which led occasionally to error (I recall you also discovered this in the test suite for Sylvester equation solvers). So, the standard test with the current Julia version is working, but the nightly test is still failing due to the hvcat issue. How long it takes to include your PR?

I observed that I have a long list of branches on my remote repository:

$ git branch -r
  dkarrasch/compathelper/new_version/2021-06-22-01-01-46-265-3581712447
  dkarrasch/compathelper/new_version/2021-06-23-00-57-02-098-2220726460
  dkarrasch/compathelper/new_version/2021-06-24-00-52-06-050-402095742
  dkarrasch/dk/types
  dkarrasch/gh-pages
  dkarrasch/massinstallaction/set-up-CI
  dkarrasch/massinstallaction/set-up-CompatHelper
  dkarrasch/massinstallaction/set-up-TagBot
  dkarrasch/master
  origin/HEAD -> origin/master
  origin/gh-pages
  origin/master

Excepting a few branches (e.g., master and remotes/origin/gh-pages), I think all the rest can be deleted. Or?

Curiously, on the GitHub interface, I see only a part of these branches (more exacly only 2 branches, no one containing dkarrasch). So, I think only you can delete those phantom branches using $ git fetch --all --prune

I would try to make a version bump. Since your modifications were substantial, I propose to make a bump to v2.1. Is this OK for you?

[dkarrasch]

Don't worry about the nightly issue. Julia triage is going to discuss whether the root cause should be reverted or the fix added. But in any case, v1.8 is not going to be released with this issue open.

Concerning the many branches, you have now probably two remotes, origin and dkarrasch. The many other branches are not mine, I don't know where all the massinstallaction stuff is coming from. I have deleted them no from my fork and locally.

As you noticed, I haven't turned the adj* fields to types in the Inverse[Generalized]SylvesterMap, because these have two and I felt it a little overkill, but one could. If not, then I think we can close this issue.

[andreasvarga]

Indeed there are two remotes:

$ git remote
dkarrasch
origin

I removed some branches with prune:

$ git remote prune  dkarrasch
Pruning dkarrasch
URL: https://github.com/dkarrasch/matrixequations.jl
 * [pruned] dkarrasch/compathelper/new_version/2021-06-22-01-01-46-265-3581712447
 * [pruned] dkarrasch/compathelper/new_version/2021-06-23-00-57-02-098-2220726460
 * [pruned] dkarrasch/compathelper/new_version/2021-06-24-00-52-06-050-402095742
 * [pruned] dkarrasch/massinstallaction/set-up-CI
 * [pruned] dkarrasch/massinstallaction/set-up-CompatHelper
 * [pruned] dkarrasch/massinstallaction/set-up-TagBot

I still don't know how to get rid of some branches, which are not mines. I believe they somehow belong to the forked repository. Here is the current status as I am seeing from my local computer:

$ git branch -r
  dkarrasch/dk/types
  dkarrasch/gh-pages
  dkarrasch/master
  origin/HEAD -> origin/master
  origin/gh-pages
  origin/master

If you want you can close this issue, to which, the above discussion is not relevant.