Comments (8)
dpo
commented on June 9, 2024
If I understand your question correctly, symmetry of the solution of your Lyapunov equation doesn't really depend on the operator, but rather on how you model the equation. Couldn't you formulate the problem by only considering tril(X) as your set of unknowns?
If I misunderstood, could you give a short example illustrating the issue you're having?
from linearoperators.jl.
andreasvarga
commented on June 9, 2024
Let consider the inverse operator inv(L) of a Lyapunov operator L(X)
=AX+XA'. This operator maps a symmetric matrix Y to a corresponding
symmetric X which satisfies the Lyapunov equation AX+XA'+Y =0. Thus
inv(L)*vec(Y) = vec(X). To determine X you need to solve the above Lyapunov
equation, and Y must be symmetric. If a specialized solver is used which
exploits symmetry, a substantial gain in efficiency can be achieved. If
such a solver underlies the product operation, then inv(L)*Y can be
computed only for a symmetric Y. However, this would lead to error when
executing Matrix(inv(L)) due to the way how the n*n columns of this matrix
are generated using columns of the identity matrix of order n*n (n is the
order of A). Excepting n columns, the rest of columns do not correspond to
symmetric matrices and the solution process would fail. Would be it
possible to extend the definition of the linear operators by specifying the
classes of structured matrices on which the operator acts?
Dominique <[email protected]> schrieb am Di., 15. Okt. 2019, 17:41:
… If I understand your question correctly, symmetry of the solution of your
Lyapunov equation doesn't really depend on the operator, but rather on how
you model the equation. Couldn't you formulate the problem by only
considering tril(X) as your set of unknowns?
If I misunderstood, could you give a short example illustrating the issue
you're having?
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
<#101?email_source=notifications&email_token=ALJDHECWGF7BYNHROM6JTXLQOXQDXA5CNFSM4JAZ5SQ2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEBJHYJI#issuecomment-542276645>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ALJDHEFGMAYWRZQQ5OHXWU3QOXQDXANCNFSM4JAZ5SQQ>
.
from linearoperators.jl.
dpo
commented on June 9, 2024
So it's the default implementation of Matrix(op::AbstractLinearOperator) that is causing you trouble. An easy solution might be to define your own operator type, say LyapunovOperator <: AbstractLinearOperator, and define the product operations accordingly. Then you can implement your own Matrix(op::LyapunovOperator) = op * Matrix(1.0I, size(op)...) (or whichever way works for you)?! Does that make sense?
from linearoperators.jl.
andreasvarga
commented on June 9, 2024
Yes, this could be an elegant solution to my problem. In the meantime, I think I found a less elegant way to manage the issue, by assuming that only the upper triangular part of X is packed in the vector x, which is then used in the product op * x. The result is the packed upper triangular part of AX+XA’. So the corresponding matrix is not n^2 times n^2 but only n(n+1)/2 times n(n+1)/2 and the special solvers can be used to implement the inverse operators.
Many thanks for your time.
Another issue: It could be helpful for your LinearOperators project a LAPACK wrapper which I implemented as an interface to the functions *lanc2.f which is instrumental to estimate efficiently the 1-norm of linear operators, and thus can serve for condition number estimation. This function is in the file lapackutil.jl in the previously mentioned package https://github.com/andreasvarga/MatrixEquations.jl
I am still working hardly towards version 1.0 to be registered, but this basic tool will not be changed (I hope) in the near future.
from linearoperators.jl.
dpo
commented on June 9, 2024
I think I found a less elegant way to manage the issue, by assuming that only the upper triangular part of X is packed in the vector x
Yes, that's what I meant in my first reply.
I'm open to adding a function for estimating the 1-norm of an operator but it cannot be based on LAPACK because we can't assume that operators are explicit matrices.
from linearoperators.jl.
andreasvarga
commented on June 9, 2024
No matrices are involved, only vectors which are the results of product
operations op*x and op'*x. So, it seems to be perfectly suited for
LinearOperators.
Dominique <[email protected]> schrieb am Mi., 16. Okt. 2019, 17:00:
… I think I found a less elegant way to manage the issue, by assuming that
only the upper triangular part of X is packed in the vector x
Yes, that's what I meant in my first reply.
I'm open to adding a function for estimating the 1-norm of an operator but
it cannot be based on LAPACK because we can't assume that operators are
explicit matrices.
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
<#101?email_source=notifications&email_token=ALJDHEEPGK55PCO2I2RN7LLQO4UBPA5CNFSM4JAZ5SQ2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEBM2ASQ#issuecomment-542744650>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ALJDHEACDNTADQ3TAPCZFXDQO4UBPANCNFSM4JAZ5SQQ>
.
from linearoperators.jl.
dpo
commented on June 9, 2024
I see. It's reverse communication. That's an idea but it would be better in the long run to rewrite it in Julia, and have all four versions (and much more) for the price of one. It's simple enough: http://www.netlib.org/lapack/explore-html/d3/d0a/dlacn2_8f_source.html
from linearoperators.jl.
andreasvarga
commented on June 9, 2024
Just an update: I formulated for my package an issue which is related to the above discussion. Simply said, I was not able to ensure the condition Matrix(T)' = Matrix(T') for an operator T which maps the upper triangular part of a symmetric matrix X to the upper triangular part of Y := AX+XA'. Your comment on this issue would be highly appreciated.
from linearoperators.jl.
Related Issues (20)
- Add a buildkite workflow
- Type inference and exponentially large compile times HOT 4
- A general QNOperator HOT 1
- Link to tutorial is broken HOT 1
- Efficient `get_nargs` HOT 3
- Any plan for including multidimensional linear operators? HOT 2
- Saddle Operator HOT 4
- Chainrules
- Separate updates for Diagonal QN operators
- Warning during the precompilation HOT 1
- Change LDLFactorizations to be a weakdep, i.e., and extension
- Remove shape function (duplicate of size) HOT 1
- Declare missing compat entries HOT 1
- Bug in derived adjoint operation? HOT 2
- `InverseLBFGSOperator` for GPU HOT 1
- `ERROR: LinearOperatorException("storage types cannot be promoted to a concrete type")` HOT 7
- Bug: adjoint(opEye()) * randn(10) HOT 2
- Review julia compatibility in readme
- ChainRulesCore error HOT 4
- Version 2.7.0 is breaking (though the version might indicate it is not) – DiagonalQN used in other packages. HOT 2
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
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.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from linearoperators.jl.