Implement SubBasis

[kalmarek]

@blegat

[codecov]

Codecov Report

Attention: Patch coverage is 73.07692% with 7 lines in your changes missing coverage. Please review.

Project coverage is 84.55%. Comparing base (bbcd649) to head (f0a91dd).

Files with missing lines	Patch %	Lines
src/bases_fixed.jl	73.07%	7 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #61      +/-   ##
==========================================
- Coverage   85.93%   84.55%   -1.39%     
==========================================
  Files          14       13       -1     
  Lines         761      777      +16     
==========================================
+ Hits          654      657       +3     
- Misses        107      120      +13     

Flag	Coverage Δ
unittests	84.55% <73.07%> (-1.39%)	⬇️

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.

[blegat]

This should be convert(I, b.parent_basis[b.supporting_elts[x]]) because we store the indices

[blegat]

This should be the list of indices of the parent_basis. For DiracBasis, since indices are the same as elements, it works as well

[kalmarek]

@blegat I was pondering this topic over and over;
I have the feeling that the whole thing should be fixed on the level of coefficients instead of basis. And maybe that's why you should provide a few test applications, so that I can see where's the crux?

At the moment what we have is

• B -- the parent basis
• A1, A2 -- subbases of B (listing indices of B) and
• algebra elements:
  • f1 contains (A1, v1) where v1 is the vector of coefficients wrt A1)
  • f2 contains (A2, v2).

We need to add/multiply f1 and f2, so we

1. construct c1 a coefficients object from v1 and A1 (i.e. the one that implements nonzero_pairs)
2. construct c2 for v2 and A2
3. pass those into our mul!/add! functions

• inside we need to construct sc::SparseCoefficients containing the result (this will be non-generic code btw. as preallocate in the generic is called only once and the result is stored directly there)

4. construct A3 -- a subbasis of B from sc and the corresponding vector v3
5. f3 needs to store now (A3, v3).

It seems that all of this can be achieved on the level of coefficients.
Why can't f1 already contain coefficients with respect to B which internally store (v1, list of corresponding indices of B)?

[blegat]

Good question. I definitely need an ExplicitBasis like SubBasis to represent a basis not necessarily used by an AlgebraElement but maybe we could consider that whenever you create an algebra element, you always use the parent ImplicitBasis. This is actually what I already do in MultivariateBases and SumOfSquares. I do all computation with implicit bases and then I just use explicit basis with coeffs.
The big blocker at the moment is the ability to represent elements from https://github.com/kalmarek/SymbolicWedderburn.jl/
We have a basis b1 of monomials of degree d and a basis b2 of monomials of degree up to 2d.
We do symmetry adapted basis for b1 -> a1 and a symmetry adapted basis for b2 -> a2.
We need to compute coeffs(a2, a1' * Q * a1).
So we should create a basis that concatenates a1 and a2 and have coeffs(a2, a1' * Q * a1) working.
If we can refactor https://github.com/kalmarek/SymbolicWedderburn.jl/blob/master/examples/sos_problem.jl using StarAlgebras then it means I can use that in SumOfSquares. At the moment, symmetry reduction in SumOfSquares has not been updated yet so that's the blocker to release a new version of SumOfSquares. Since in SumOfSquares, everything need to implement the StarAlgebras interface now, I need to implement coeffs(a2, a1' * Q * a1) using StarAlgebras so once we have a StarAlgebras basis for SymbolicWedderburn all lights will be green.
So the crux is https://github.com/kalmarek/SymbolicWedderburn.jl/blob/master/examples/sos_problem.jl / kalmarek/SymbolicWedderburn.jl#78 !

[kalmarek]

@blegat, I'm not sure what do you mean by

If we can refactor https://github.com/kalmarek/SymbolicWedderburn.jl/blob/master/examples/sos_problem.jl using StarAlgebras then it means I can use that in SumOfSquares.

Nowhere there do we form a product coeffs(a2, a1'*Q*a1) and that code after small tweaks should work with the newest StarAlgebras.

Here's some newly updated code that make use of DiracBasis and switches to FixedBasis only when particular (fixed) dimension relaxation is required:

import StarAlgebras as SA
let X = .... # AlgebraElement for a StarAlgebra of FPMonoid (infinite dimensional, DiracBasis)
    A = parent(elts)
    unit = one(A)

    # compute somehow a set of elements from `basis(A)` where dual sos-problem for X makes sense.
    linsupp, sizes = linear_supp(X, unit)
    B = SA.FixedBasis(linsupp, SA.mstructure(SA.basis(A)))
    
    # length(word(x)) corresponds to degree of polynomials
    N = maximum(x -> sizes[length(FPMonoids.word(x))÷2], linsupp)
        
    mt = let psd_basis = @view linsupp[1:N] # linsupp is ordered by word-length
        # here we use DiracMStructure of `B` and `mt` is `Matrix{<:Integer}`
        mt = [B[star(x)*y] for x in psd_basis, y in psd_basis]
    end

    model = JuMP.Model()
    # y is a vector of variables in basis `B`!
    y = @variable model y[1:length(B)]
    @objective model Max dot(coeffs(X, B), y) # max dot(X, y)
    @constraint model dot(coeffs(unit, B), y) == 1 # given a normalized `y`
    let P = @view y[mt]
        @constraint model P in PSDCone() # and mt defines a psd matrix
    end

    return model
end

see also https://github.com/kalmarek/QuantumStuff.jl/blob/mk/new_fpmonoids/src/psd_problems.jl

[blegat]

Thanks, I'll play around with while trying to update jump-dev/SumOfSquares.jl#375

[blegat]

Meeting notes:
q is SOS
∃ Q PSD s.t.
q == g -> on some basis zero_basis
where g := gram_basis' * Q * gram_basis
coeffs(zero_basis, gram_basis' * Q * gram_basis)

zero_basis -> Any SA.ExplicitBasis
gram_basis -> Any SA.ExplicitBasis
For all a, b in gram_basis, a*b ∈ zero_basis

coeffs(zero_basis, q) == coeffs(zero_basis, gram_basis' * Q * gram_basis)

How to compute coeffs(zero_basis, gram_basis' * Q * gram_basis) ???

function quadratic_form(Q, basis)
    p = zero(implicit_basis(zero_basis))
    MA.operate!(..., p, gram_basis' * Q * gram_basis)
    return p
end
coeffs(zero_basis, quadratic_form(gram_matrix.Q, gram_matrix.basis))

Symmetry case:
gram_basis Fixed basis containing polynomials or more generally AlgebraElement{FullBasis{B}}
B could be Monomial or B could be Chebyshev

MA.operate!(..., p, gram_basis' * Q * gram_basis)
for row, col
   MA.operate!(...,
       p::AlgebraElement{FullBasis{B}},
       1 * star(gram_basis[row])::SymmetryAdapted{B},
       Q[row, col] * gram_basis[col]::SymmetryAdapted{B},
   )
end
MA.operate!(..., p::AlgebraElement{FullBasis{B}}, ::AlgebraElement{FullBasis{B}})

from SOS

function MA.operate!(
    op::SA.UnsafeAddMul{typeof(*)},
    p::SA.AlgebraElement,
    g::GramMatrix,
    args::Vararg{Any,N},
) where {N}
    for row in eachindex(g.basis)
        row_star = SA.star(g.basis[row])
        for col in eachindex(g.basis)
            MA.operate!(
                op,
                p,
                MB.term_element(g.Q[row, col], identity(algebra(p))),
                MB.convert_basis(row_star, SA.algebra(p)),
                MB.convert_basis(g.basis[col], SA.algebra(p)),
                args...,
            )
        end
    end
    return p
end

Conclusion, we should have something like this in StarAlgebras: operate!(::UnsafeAddMul, ::AlgebraElement, QuadraticForm).

[kalmarek]

struct QuadraticForm{T}
    Q::T
end

basis(Q::QuadraticForm) = basis(Q.Q) # ExplicitBasis
Base.getindex(Q, i::T, j::T) where {T} = Q[basis(Q)[i], basis(Q)[j]]

function Base.getindex(Q, i::Integer, j::Integer)
    # returns the q_ij
    return Q[i, j]
end

function MA.operate!(
    op::UnsafeAddMul,
    p::T, # implicit basis
    Q::QuadraticForm,
)
    for b1 in basis(Q)
        b1′ = star(b1)
        for b2 in basis(Q)
            MA.operate!(op, p, (1, b1′), (Q[i, j], b2))
        end
    end
    return p
end