Implement quadratic forms

[kalmarek]

@blegat I hacked the form adding cfs kwarg to operate! on UnsafeAddMul;

QuadraticForm(Q) is a simple quadratic form.
QuadraticForm{star}(Q) is a sesqui-linear form (but if star === identity on your basis it will not perform additional starring).

---

A few gripes:

1. I don't like all of those args as they bring little value and only obfuscate the code
2. debugging call stack is a nightmare for anyone not intimately familiar with MA. It's a technical debt that will bite us in our arses;
3. operate!(op::UnsafeAddMul, res, c) should become operate!(op::UnsafeAdd, res, c)
4. we're gaining nothing by passing around the args (only a recursive call)

I propose implementing it via loop (of length known to the compiler) via UnsafeLAMul! (or whatever name) with mul!(C, A, B, α, β) performing in-place A·B·α + C·β.
As we don't currently need β, I'd vote for 4-arg mul!(C, A, B, α) until β is needed.

Advantages:

1. No recursive calls
2. Easy to un-hack cfs as α
3. vararg can be implemented with a simple loop on top;

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 83.56%. Comparing base (bbcd649) to head (51d67fe).

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #62      +/-   ##
==========================================
- Coverage   85.93%   83.56%   -2.38%     
==========================================
  Files          14       14              
  Lines         761      785      +24     
==========================================
+ Hits          654      656       +2     
- Misses        107      129      +22     

Flag	Coverage Δ
unittests	83.56% <100.00%> (-2.38%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[kalmarek]

@blegat there are some problems with allocations of simple
@allocated(MA.operate_to!(d, *, 2, a)) == 0 ??

What broke it? MA allocation stuff seems really fragile...

[blegat]

I do need the args... in case the user uses the domain argument in SumOfSquares, it's the Putinar certificate.

[blegat]

Suggested change

	struct Gramm{T,B}
	struct Gram{T,B}

[blegat]

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

[kalmarek]

• b1★ and b2 come from basis(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 as b1 and b2, it's enough that their keys are in the basis of res (e.g. they are defined over some fixed basis, but res 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

[blegat]

b1★ and b2 come from basis(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 as b1 and b2, it's enough that their keys are in the basis of res (e.g. they are defined over some fixed basis, but res 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 ;)

[kalmarek]

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

[blegat]

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! :)

[kalmarek]

have a look now;

MA would really benefit from proper documentation (some design principles, maybe a tutorial, show-casting the call stack etc.)

[blegat]

I propose implementing it via loop (of length known to the compiler) via UnsafeLAMul! (or whatever name) with mul!(C, A, B, α, β) performing in-place A·B·α + C·β.
As we don't currently need β, I'd vote for 4-arg mul!(C, A, B, α) until β is needed.

I'm ok with not using MA for UnsafeAddMul, until there is a canonical! notion in MA, we don't gain much by using the MA interface for it.
mul! will just overwrite and erase C while here we need to add to C, leaving the terms that are already in C. It is also unsafe because there are duplicates so it should be unsafe_add_mul!.
I need at least 3 algebra elements for the putinar certificate so we could have unsafe_add_mul!(*, C, A, B, args...)! and unsafe_add_mul_with_constant!(*, C, α, A, B, args...).

[blegat]

I actually already tried adding scalars in the mix in
#34 (comment)
But I agree it's best to stay with at most one scalar and have an explicit place for it (keyword argument or second argument).

[kalmarek]

I propose implementing it via loop (of length known to the compiler) via UnsafeLAMul! (or whatever name) with mul!(C, A, B, α, β) performing in-place A·B·α + C·β.
As we don't currently need β, I'd vote for 4-arg mul!(C, A, B, α) until β is needed.

I'm ok with not using MA for UnsafeAddMul, until there is a canonical! notion in MA, we don't gain much by using the MA interface for it. mul! will just overwrite and erase C while here we need to add to C, leaving the terms that are already in C. It is also unsafe because there are duplicates so it should be unsafe_add_mul!. I need at least 3 algebra elements for the putinar certificate so we could have unsafe_add_mul!(*, C, A, B, args...)! and unsafe_add_mul_with_constant!(*, C, α, A, B, args...).

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 unsafe_add_scalar_mul_vararg!?

[blegat]

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 unsafe_add_scalar_mul_vararg!?

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 Number.

[blegat]

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)

[blegat]

I'm going to make a PR on SumOfSquares using this to double check

[blegat]

Maybe drop Julia v1.6 and add a v1.10 as minimum requirement, it's the new LTS anyway

[blegat]

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 😇

[kalmarek]

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

[kalmarek]

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

[blegat]

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

[blegat]

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))

[kalmarek]

p = zero(...)
MA.operate!(UnsafeAddMul(...), p, SA.QuadraticForm(q))

what is multiplied here? do you mean p += AlgebraElement(qf, parent(p)?

p = zero(...)
for (q, g) in zip(...)
    MA.operate!(UnsafeAddMul(...), p, SA.QuadraticForm(q), g)
end

I think doing this will be a huge performance bottleneck.
you will end up in nterms(q)^2*nterms(g) terms in p.
How many terms will there be after canonicalization of qf?

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) nterms(q)^2 is ~2e7, while nterms(qf) is <1e7 (after canonical). Do you really want to pay double (times nterms(g)) the price?