-
Notifications
You must be signed in to change notification settings - Fork 3
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implement quadratic forms #62
base: main
Are you sure you want to change the base?
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #62 +/- ##
==========================================
- Coverage 85.93% 83.31% -2.63%
==========================================
Files 14 14
Lines 761 791 +30
==========================================
+ Hits 654 659 +5
- Misses 107 132 +25
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. |
@blegat there are some problems with allocations of simple What broke it? MA allocation stuff seems really fragile... |
I do need the |
src/mstructures.jl
Outdated
for (i, b1) in pairs(basis(Q)) | ||
b1★ = ε(b1) | ||
for (j, b2) in pairs(basis(Q)) | ||
MA.operate!(op, res, coeffs(b1★), coeffs(b2); cfs = Q[i, j]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why is res
already a coefficient ? I would simply call MA.operate!(op, res, b1★, b2; cfs = Q[i, j])
and let it follow the normal call stack that will then unwrap the algebra elements calling coeffs
but first checking they are of the same algebra which is safer
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
b1★
andb2
come frombasis(Q)
, i.e. the assumption is already fulfilled- to go through the regular call stack we'd need to support
cfs
kwarg throughout the whole stack; maybe it is sensible, or maybe it isn't; res
doesn't have to be in the same algebra asb1
andb2
, it's enough that theirkeys
are in the basis ofres
(e.g. they are defined over some fixed basis, butres
is done over implicit one).
But the actual true reason is that it's already too complicated for me to follow the call stack mentally ;D There are too many operate!
, operate_to!
methods which I don't even know what's the difference and where are they located in the stack. So I took shortcut just to make this work
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
b1★
andb2
come frombasis(Q)
, i.e. the assumption is already fulfilled
How do we know it is ?
to go through the regular call stack we'd need to support
cfs
kwarg throughout the whole stack; maybe it is sensible, or maybe it isn't;
Not the whole call stack, only the unsafe add mul one. It should be supported for unsafe add mul with algebra elements and with coefficients.
res
doesn't have to be in the same algebra asb1
andb2
, it's enough that theirkeys
are in the basis ofres
(e.g. they are defined over some fixed basis, butres
is done over implicit one).
If we have an unsafe add mul call that writes to an implicit, it could indeed support that some of the arguments are in explicit form indeed
But the actual true reason is that it's already too complicated for me to follow the call stack mentally ;D There are too many
operate!
,operate_to!
methods which I don't even know what's the difference and where are they located in the stack. So I took shortcut just to make this work
We could not use the MA API for the unsafe add mul method but for the other one, it's important we still use it because we gain all the optimization that were done with the rewriting for JuMP at the cost of the added complexity ;)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
basis here is of a star algebra (considered as a vector space), hence spanned by vectors indexed by elements of the algebra. But this is only a soft assumption, what we need to know is that the multiplicative structure of the basis of res
handles the keys of products of basis(Q)
well.
Not the whole call stack, only the unsafe add mul one. It should be supported for unsafe add mul with algebra elements and with coefficients.
I went with a 4-argument unsafe add mul with the form MA.opearte_to!(res, ::UnsafeAddMul, A, B, α)
which computes α·A·B
storing the result in res
. Have a look at the changes.
Since we always pretend that res are truly mutable I changed all the methods to MA.operate_to!(res, ...)
. I still don't understand the difference of operate!(::UnsafeMulAdd, res, ...)
and operate_to!(res, ::UnsafeMulAdd, ...)
.
Well the benefit right now is that I'm afraid of touching those MA.something
methods. They feel to me like typesetting latex table :P
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's not operate_to!
because operate_to!
overwrites the result while we just want to append to it (hence the Add
in UnsafeAddMul
. So it's really operate!
:)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
have a look now;
MA would really benefit from proper documentation (some design principles, maybe a tutorial, show-casting the call stack etc.)
I'm ok with not using |
I actually already tried adding scalars in the mix in |
This starts looking very much like C, where it's the best to encode the signature in the name of the function... did you mean |
And the beauty is, you can do C-like code without having to leave the language ;) We can have the same function with and without the constant if we have a convention that the constant appears before the factors but then we need to have a convention on what is the type of the constant or what is the type of the coefficients. Here, the constant could be a JuMP expression so we cannot assume that it's a |
In the meantime, I'm making progress on jump-dev/SumOfSquares.jl#375 to clarify what the target is (we consume the rope by burning it on both sides) |
I'm going to make a PR on SumOfSquares using this to double check |
Maybe drop Julia v1.6 and add a v1.10 as minimum requirement, it's the new LTS anyway |
basis(qf::QuadraticForm) = basis(qf.Q) | ||
Base.getindex(qf::QuadraticForm, i::T, j::T) where {T} = qf.Q[i, j] | ||
|
||
function MA.operate_to!( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This should be operate!(ms, res, Q, args...)
. For Putinar, I need to add several quadratic form so it won't work if you zero
. And I need the extra args 😇
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do I read correctly, that you need to add several quadratic forms together?
E.g.
function MA.operate_to!(res, ms::MultiplicativeStructure, Qs::Vararg{<:QuadraticForm})
MA.operate!(zero, res)
for Q in Qs
MA.operate!(ms, res, Q)
end
MA.operate!(canonical, res)
return res
end
function MA.operate!(ms::MultiplicativeStructure, res, Q::QuadraticForm{T,ε}) where {T,ε}
op = UnsafeAddMul(ms)
for (i, b1) in pairs(basis(Q))
b1★ = ε(b1)
for (j, b2) in pairs(basis(Q))
MA.operate!(op, res, coeffs(b1★), coeffs(b2), Q[i, j])
end
end
return res
end
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It would be much easier to push this over the finish line had you provided the use-cases you have in mind for SoS as I asked many a times... :D
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm as lost in the design space as you are ^^
MA.operate_to!(res, ms::MultiplicativeStructure, Qs::Vararg{<:QuadraticForm})
doesn't work because it's a sum, MA.operate_to!(res, ::typeof(+), Qs::Vararg{<:QuadraticForm})
is correct but I won't use it because I will multiply each graph by a poynomial to do putinar
So the use case is
p = zero(...)
for (q, g) in zip(...)
MA.operate!(UnsafeAddMul(...), p, QuadraticForm(q), g)
end
So the use case is p = zero(...)
for (q, g) in zip(...)
MA.operate!(UnsafeAddMul(...), p, SA.QuadraticForm(q), g)
end The other use case when it's not putinar is: p = zero(...)
MA.operate!(UnsafeAddMul(...), p, SA.QuadraticForm(q)) |
what is multiplied here? do you mean
I think doing this will be a huge performance bottleneck. I honestly doubt it will be beneficial for large scale examples. In one of my cases (which is worse by a lot than any commutative case) |
@blegat I added a 3-argument A = SA.StarAlgebra(...)
Q = SA.QuadraticForm(...)
res = zero(eltype(Q), A)
p = sum(basis(Q)) # a 'random' element from A
MA.operate_to!(coeffs(res), SA.mstructure(A), Q, p)
@test res == A(Q) * p
|
Yes
It may be, the user can canonicalize after adding each term if needed. We could also have a
Just tested, it almost works. I have the issue that I do |
src/quadratic_form.jl
Outdated
for (i, b1) in pairs(basis(Q)) | ||
b1★ = ε(b1) | ||
for (j, b2) in pairs(basis(Q)) | ||
MA.operate!(op, res, coeffs(b1★), coeffs(b2), Q[i, j]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
So every basis element should implement coeffs
? It's weird because it's not really an AlgebraElement
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What if the elements of the basis are just group elements
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If we see AlgebraElement, this implementation is correct. Otherwise we should multiply b1 and b2 using the multiplicative structure
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It can just return coeffs(x) = (x, true)
if you like;
They are not just group elements, they are Diracs
at group elements. The proper formalism is to talk about functions (on groups or other objects) instead of formal sums.
So far it seems that you can place in a basis something sortable (unless you canonicalize!
yourself) that implements coeffs
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you use arbitrary SA.ExplicitBasis
it won't work right ? For example you place group elements in the basis and coeffs
is not part of GroupsCore
API.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would just call some new internal function and distinguish two cases with multiple dispatch
- Two
AlgebraElement
-> then do what is currently done - Otherwise, just multiply them with the multiplicative structure
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You don't just place group elements in basis, you need to pick a multiplicative structure. It could be e.g. DiracMStructure
and then Diracs
are in your basis (they aren't AlgebraElements
either).
Does your element behave like coefficients? then we could define SA.coeffs(x::YourTypeWhichBehavesLikeAbstractCoeffs) = x
to push x
down the call?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
but in general, I'm not convinced that this quadratic form is a good design. Let's brainstorm for a bit where should all coeffs
calls go.
Can you formulate a concrete example specifying:
- What will be in the basis
- what is the mstructure
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok let me formulate a use case in the main thread
Co-authored-by: Benoît Legat <[email protected]>
What I meant is with the current code the user cannot escape creating all of those redundant terms even if the storage is canonicalizing on the spot (via dictionary or something). I'm not convinced that the triple product (Putinar) is even beneficial to have. Small example will run fast either with or without allocations; large examples will suffer from algorithmic complexity. |
I see, you mean that we should first transform the gram matrix into an |
@blegat What is left here?
|
) | ||
# ABC::Vararg{AlgebraElement,N}, | ||
# ) where {N} | ||
MA.operate!(mul, coeffs(res), coeffs(A), coeffs(B), α) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why don't we check that they have the same parent here ?
@@ -115,7 +115,7 @@ function LinearAlgebra.dot(w::AbstractVector, ac::AbstractCoefficients) | |||
if isempty(values(ac)) | |||
return zero(MA.promote_sum_mul(eltype(w), value_type(ac))) | |||
else | |||
return sum(w[i] * star(v) for (i, v) in nonzero_pairs(ac)) | |||
return sum(w[i] * v for (i, v) in nonzero_pairs(ac)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why don't we take the conjugation?
end | ||
MA.operate!(zero, res) | ||
MA.operate!(UnsafeAddMul(ms), res, args...) | ||
res = MA.operate!(UnsafeAddMul(ms), res, args...) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
With operate!
you are guaranteed that the argument res
is the same as what is returned so no need to get the result
end | ||
|
||
function MA.operate!( | ||
op::UnsafeAddMul, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Now that we have UnsafeAdd
, we should use UnsafeAdd
here
|
||
function MA.operate_to!( | ||
res, | ||
ms::MultiplicativeStructure, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What's the plan for this ms
? Is the plan to allow res
to have a different basis compared to Q and ms is the rule to do the mapping ? If yes, this should be commented on. Otherwise, just use copy
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
yeah, this signature needs to change. What we want here is to say:
Take this quadratic form (which comes with its own basis) and realize it algebraically storing the result in res
.
Is copy
the correct convention? store
of copyto
could also work;
how about MA.operate_to!(res, Q)
where the quadratic form itself "operates"?
(this ms
is an artifact of the past.)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There is precedent for copy
in
https://github.com/jump-dev/MutableArithmetics.jl/blob/master/src/implementations/BigInt.jl#L20
It's weird that it's not the same type though so it might not be the right one.
Having Q
be the operator is weird also. We could use UnsafeAdd
here, it's like the unary +
.
|
||
function MA.operate_to!( | ||
res, | ||
ms::MultiplicativeStructure, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ms::MultiplicativeStructure, | |
::typeof(copy), |
Q::QuadraticForm, | ||
) | ||
MA.operate!(zero, res) | ||
MA.operate!(UnsafeAddMul(ms), res, Q) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
MA.operate!(UnsafeAddMul(ms), res, Q) | |
MA.operate!(UnsafeAdd(), res, Q) |
Let's focus on making this work, this is what is needed for SumOfSquares. b1 = # explicit basis containing monomials or group elements
A = # algebra over the implicit version of b1
b2 = # explicit basis containing `AlgebraElement`s over A
Q1 = # QuadraticForm with b1
Q2 = # QuadraticForm with b2
AlgebraElement(Q1, A) # should be over A
AlgebraElement(Q2, A) # should be over A |
@blegat I hacked the form adding
cfs
kwarg tooperate!
onUnsafeAddMul
;QuadraticForm(Q)
is a simple quadratic form.QuadraticForm{star}(Q)
is a sesqui-linear form (but ifstar === identity
on your basis it will not perform additional starring).A few gripes:
args
as they bring little value and only obfuscate the codeoperate!(op::UnsafeAddMul, res, c)
should becomeoperate!(op::UnsafeAdd, res, c)
I propose implementing it via loop (of length known to the compiler) via
UnsafeLAMul!
(or whatever name) withmul!(C, A, B, α, β)
performing in-placeA·B·α + C·β
.As we don't currently need
β
, I'd vote for 4-argmul!(C, A, B, α)
untilβ
is needed.Advantages:
cfs
asα