From 2b8f0fb7fa8c4ae8f5ef7a6c91881251b2c04483 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 15 Dec 2024 15:14:10 +0100 Subject: [PATCH] Fix symmetry reduction (#375) * Fix symmetry reduction * Fix * Fixes * Fixes * Fixes * Reenable Symmetry * Fixes * Fix * Remove commented code * Fix format * Updates * Update * Fix * Fix format --- docs/make.jl | 2 +- docs/src/tutorials/Symmetry/cyclic.jl | 65 ++++++----- docs/src/tutorials/Symmetry/dihedral.jl | 74 +++++++------ docs/src/tutorials/Symmetry/even_reduction.jl | 11 +- .../Symmetry/permutation_symmetry.jl | 47 ++++---- src/Certificate/Symmetry/Symmetry.jl | 1 + src/Certificate/Symmetry/wedderburn.jl | 101 +++++++++++++----- src/Certificate/ideal.jl | 2 +- src/gram_matrix.jl | 63 +++++++++-- 9 files changed, 241 insertions(+), 125 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index a6073507d..f56db8f8e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,7 +13,7 @@ const _TUTORIAL_SUBDIR = [ "Other Applications", "Noncommutative and Hermitian", "Sparsity", - #"Symmetry", + "Symmetry", "Extension", ] diff --git a/docs/src/tutorials/Symmetry/cyclic.jl b/docs/src/tutorials/Symmetry/cyclic.jl index f7088f4b0..77c4280c9 100644 --- a/docs/src/tutorials/Symmetry/cyclic.jl +++ b/docs/src/tutorials/Symmetry/cyclic.jl @@ -51,6 +51,7 @@ end # `x_1^3*x_2*x_3^4` into `x_1^4*x_2^3*x_3`. import MultivariatePolynomials as MP +import MultivariateBases as MB using SumOfSquares @@ -83,12 +84,24 @@ Symmetry.SymbolicWedderburn.action(action, g, poly) # Let's now find the minimum of `p` by exploiting this symmetry. -import CSDP -solver = CSDP.Optimizer +import Clarabel +solver = Clarabel.Optimizer model = Model(solver) @variable(model, t) @objective(model, Max, t) pattern = Symmetry.Pattern(G, action) +basis = MB.explicit_basis(MB.algebra_element(poly - t)) +using SymbolicWedderburn +summands = SymbolicWedderburn.symmetry_adapted_basis( + Float64, + pattern.group, + pattern.action, + basis, + semisimple = true, +) + +gram_basis = SumOfSquares.Certificate.Symmetry._gram_basis(pattern, basis, Float64) + con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern) optimize!(model) @test termination_status(model) == MOI.OPTIMAL #src @@ -98,21 +111,19 @@ solution_summary(model) # Let's look at the symmetry adapted basis used. gram = gram_matrix(con_ref).blocks #src -@test length(gram) == 3 #src +@test length(gram) == 2 #src @test gram[1].Q ≈ [0 0; 0 2] #src -@test length(gram[1].basis.polynomials) == 2 #src -@test gram[1].basis.polynomials[1] == 1 #src -@test gram[1].basis.polynomials[2] ≈ -sum(x)/√3 #src +polys = gram[1].basis.bases[].elements #src +@test length(polys) == 2 #src +@test polys[1] ≈ 1 #src +@test polys[2] ≈ -sum(x)/√3 #src @test gram[2].Q ≈ [0.5;;] #src -@test length(gram[2].basis.polynomials) == 1 #src -@test gram[2].basis.polynomials[1] ≈ (x[1] + x[2] - 2x[3])/√6 #src -@test gram[3].Q == gram[2].Q #src -@test length(gram[3].basis.polynomials) == 1 #src -@test gram[3].basis.polynomials[1] ≈ (x[1] - x[2])/√2 #src -for gram in gram_matrix(con_ref).blocks - println(gram.basis.polynomials) - display(gram.Q) -end +@test length(gram[2].basis.bases) == 2 #src +polys = gram[2].basis.bases[1].elements #src +@test polys[] ≈ (x[1] + x[2] - 2x[3])/√6 #src +polys = gram[2].basis.bases[2].elements #src +@test polys[] ≈ (x[1] - x[2])/√2 #src +gram_matrix(con_ref) # Let's look into more details at the last two elements of the basis. @@ -137,8 +148,6 @@ b = √3/2 # In fact, these last two basis comes from the real decomposition of a complex one. -import CSDP -solver = CSDP.Optimizer model = Model(solver) @variable(model, t) @objective(model, Max, t) @@ -153,19 +162,19 @@ solution_summary(model) gram = gram_matrix(con_ref).blocks #src @test length(gram) == 3 #src @test gram[1].Q ≈ [0 0; 0 2] #src -@test length(gram[1].basis.polynomials) == 2 #src -@test gram[1].basis.polynomials[1] == 1 #src -@test gram[1].basis.polynomials[2] ≈ -sum(x)/√3 #src +polys = gram[1].basis.bases[].elements #src +@test length(polys) == 2 #src +@test polys[1] ≈ 1 #src +@test polys[2] ≈ -sum(x)/√3 #src @test gram[2].Q ≈ [0.5;;] rtol = 1e-6 #src -@test length(gram[2].basis.polynomials) == 1 #src -@test gram[2].basis.polynomials[1] ≈ (basis[1] - basis[2] * im) / √2 #src +polys = gram[2].basis.bases[].elements #src +@test length(polys) == 1 #src +@test polys[] ≈ (basis[1] + basis[2] * im) / √2 #src @test gram[3].Q ≈ [0.5;;] rtol = 1e-6 #src -@test length(gram[3].basis.polynomials) == 1 #src -@test gram[3].basis.polynomials[1] ≈ (basis[1] + basis[2] * im) / √2 #src -for gram in gram_matrix(con_ref).blocks - println(gram.basis.polynomials) - display(gram.Q) -end +polys = gram[3].basis.bases[].elements #src +@test length(polys) == 1 #src +@test polys[] ≈ (basis[1] - basis[2] * im) / √2 #src +gram_matrix(con_ref) # We can see that the real invariant subspace was in fact coming from two complex conjugate complex invariant subspaces: diff --git a/docs/src/tutorials/Symmetry/dihedral.jl b/docs/src/tutorials/Symmetry/dihedral.jl index 7802a8b65..94e85867e 100644 --- a/docs/src/tutorials/Symmetry/dihedral.jl +++ b/docs/src/tutorials/Symmetry/dihedral.jl @@ -8,7 +8,6 @@ # *Symmetry groups, semidefinite programs, and sums of squares*. # Journal of Pure and Applied Algebra 192.1-3 (2004): 95-128. - using Test #src import Random #src Random.seed!(0) #src @@ -134,9 +133,9 @@ end # We can exploit this symmetry for reducing the problem using the `SymmetricIdeal` certificate as follows: -import CSDP +import Clarabel function solve(G) - solver = CSDP.Optimizer + solver = Clarabel.Optimizer model = Model(solver) @variable(model, t) @objective(model, Max, t) @@ -148,37 +147,44 @@ function solve(G) g = gram_matrix(con_ref).blocks #src - @test length(g) == 5 #src - @test g[1].basis.polynomials == [y^3, x^2*y, y] #src - @test g[2].basis.polynomials == [-x^3, -x*y^2, -x] #src - for i in 1:2 #src - I = 3:-1:1 #src - Q = g[i].Q[I, I] #src - @test size(Q) == (3, 3) #src - @test Q[2, 2] ≈ 1 rtol=1e-2 #src - @test Q[1, 2] ≈ 5/8 rtol=1e-2 #src - @test Q[2, 3] ≈ -1 rtol=1e-2 #src - @test Q[1, 1] ≈ 25/64 rtol=1e-2 #src - @test Q[1, 3] ≈ -5/8 rtol=1e-2 #src - @test Q[3, 3] ≈ 1 rtol=1e-2 #src - end #src - @test length(g[3].basis.polynomials) == 2 #src - @test g[3].basis.polynomials[1] == 1.0 #src - @test g[3].basis.polynomials[2] ≈ -(√2/2)x^2 - (√2/2)y^2 #src - @test size(g[3].Q) == (2, 2) #src - @test g[3].Q[1, 1] ≈ 7921/4096 rtol=1e-2 #src - @test g[3].Q[1, 2] ≈ 0.983 rtol=1e-2 #src - @test g[3].Q[2, 2] ≈ 1/2 rtol=1e-2 #src - @test g[4].basis.polynomials == [x * y] #src - @test size(g[4].Q) == (1, 1) #src - @test g[4].Q[1, 1] ≈ 0 atol=1e-2 #src - @test length(g[5].basis.polynomials) == 1 #src - @test g[5].basis.polynomials[1] ≈ (√2/2)x^2 - (√2/2)y^2 #src - @test size(g[5].Q) == (1, 1) #src - @test g[5].Q[1, 1] ≈ 0 atol=1e-2 #src - for g in gram_matrix(con_ref).blocks - println(g.basis.polynomials) - end + @test length(g) == 4 #src + @test length(g[4].basis.bases) == 2 #src + polys = g[4].basis.bases[1].elements #src + @test length(polys) == 3 #src + @test polys[1] ≈ y^3 #src + @test polys[2] ≈ x^2*y #src + @test polys[3] ≈ y #src + polys = g[4].basis.bases[2].elements #src + @test length(polys) == 3 #src + @test polys[1] ≈ -x^3 #src + @test polys[2] ≈ -x*y^2 #src + @test polys[3] ≈ -x #src + I = 3:-1:1 #src + Q = g[4].Q[I, I] #src + @test size(Q) == (3, 3) #src + @test Q[2, 2] ≈ 1 rtol=1e-2 #src + @test Q[1, 2] ≈ 5/8 rtol=1e-2 #src + @test Q[2, 3] ≈ -1 rtol=1e-2 #src + @test Q[1, 1] ≈ 25/64 rtol=1e-2 #src + @test Q[1, 3] ≈ -5/8 rtol=1e-2 #src + @test Q[3, 3] ≈ 1 rtol=1e-2 #src + polys = g[1].basis.bases[].elements #src + @test length(polys) == 2 #src + @test polys[1] ≈ 1.0 #src + @test polys[2] ≈ -(√2/2)x^2 - (√2/2)y^2 #src + @test size(g[1].Q) == (2, 2) #src + @test g[1].Q[1, 1] ≈ 7921/4096 rtol=1e-2 #src + @test g[1].Q[1, 2] ≈ 0.983 rtol=1e-2 #src + @test g[1].Q[2, 2] ≈ 1/2 rtol=1e-2 #src + polys = g[2].basis.bases[].elements #src + @test polys[] ≈ x * y #src + @test size(g[2].Q) == (1, 1) #src + @test g[2].Q[1, 1] ≈ 0 atol=1e-2 #src + polys = g[3].basis.bases[].elements #src + @test polys[] ≈ (√2/2)x^2 - (√2/2)y^2 #src + @test size(g[3].Q) == (1, 1) #src + @test g[3].Q[1, 1] ≈ 0 atol=1e-2 #src + gram_matrix(con_ref) end solve(G) diff --git a/docs/src/tutorials/Symmetry/even_reduction.jl b/docs/src/tutorials/Symmetry/even_reduction.jl index 22a3b2b37..856878967 100644 --- a/docs/src/tutorials/Symmetry/even_reduction.jl +++ b/docs/src/tutorials/Symmetry/even_reduction.jl @@ -31,8 +31,8 @@ G = PermGroup([perm"(1,2)"]) # We can exploit the symmetry as follows: -import CSDP -solver = CSDP.Optimizer +import Clarabel +solver = Clarabel.Optimizer model = Model(solver) @variable(model, t) @objective(model, Max, t) @@ -45,6 +45,9 @@ value(t) # We indeed find `-1`, let's verify that symmetry was exploited: @test length(gram_matrix(con_ref).blocks) == 2 #src -@test gram_matrix(con_ref).blocks[1].basis.polynomials == [1, x^2] #src -@test gram_matrix(con_ref).blocks[2].basis.polynomials == [x] #src +polys = gram_matrix(con_ref).blocks[1].basis.bases[].elements #src +@test polys[1] ≈ 1 #src +@test polys[2] ≈ x^2 #src +polys = gram_matrix(con_ref).blocks[2].basis.bases[].elements #src +@test polys[] ≈ x #src gram_matrix(con_ref) diff --git a/docs/src/tutorials/Symmetry/permutation_symmetry.jl b/docs/src/tutorials/Symmetry/permutation_symmetry.jl index a24ecec24..f2dd41f9e 100644 --- a/docs/src/tutorials/Symmetry/permutation_symmetry.jl +++ b/docs/src/tutorials/Symmetry/permutation_symmetry.jl @@ -32,8 +32,8 @@ G = PermGroup([perm"(1,2,3,4)"]) # We can use this certificate as follows: -import CSDP -solver = CSDP.Optimizer +import Clarabel +solver = Clarabel.Optimizer model = Model(solver) @variable(model, t) @objective(model, Max, t) @@ -45,25 +45,28 @@ value(t) # We indeed find `-1`, let's verify that symmetry was exploited: -g = gram_matrix(con_ref).blocks #src -@test length(g) == 4 #src -@test length(g[1].basis.polynomials) == 2 #src -@test g[1].basis.polynomials[1] == 1.0 #src -@test g[1].basis.polynomials[2] ≈ -0.5 * sum(x) #src -@test size(g[1].Q) == (2, 2) #src -@test g[1].Q[1, 1] ≈ 1.0 atol=1e-6 #src -@test g[1].Q[1, 2] ≈ -1.0 atol=1e-6 #src -@test g[1].Q[2, 2] ≈ 1.0 atol=1e-6 #src -@test length(g[2].basis.polynomials) == 1 #src -@test g[2].basis.polynomials[1] ≈ (x[2] - x[4]) / √2 #src +gram = gram_matrix(con_ref).blocks #src +@test length(gram) == 3 #src +polys = gram[1].basis.bases[].elements #src +@test length(polys) == 2 #src +@test polys[1] ≈ 1 #src +@test polys[2] ≈ -0.5 * sum(x) #src +@test size(gram[1].Q) == (2, 2) #src +@test gram[1].Q[1, 1] ≈ 1.0 atol=1e-6 #src +@test gram[1].Q[1, 2] ≈ -1.0 atol=1e-6 #src +@test gram[1].Q[2, 2] ≈ 1.0 atol=1e-6 #src +@test length(gram[2].basis.bases) == 2 #src +polys = gram[2].basis.bases[1].elements #src +@test length(polys) == 1 #src +@test polys[] ≈ (x[2] - x[4]) / √2 #src @test size(g[2].Q) == (1, 1) #src @test g[2].Q[1, 1] ≈ 1.0 atol=1e-6 #src -@test length(g[3].basis.polynomials) == 1 #src -@test g[3].basis.polynomials[1] ≈ (x[1] - x[3]) / √2 #src -@test size(g[3].Q) == (1, 1) #src -@test g[3].Q[1, 1] ≈ 1.0 atol=1e-6 #src -@test length(g[4].basis.polynomials) == 1 #src -@test g[4].basis.polynomials[1] ≈ (x[1] - x[2] + x[3] - x[4]) / 2 #src -@test size(g[4].Q) == (1, 1) #src -@test g[4].Q[1, 1] ≈ 1.0 atol=1e-6 #src -gram_matrix(con_ref).blocks +polys = gram[2].basis.bases[2].elements #src +@test length(polys) == 1 #src +@test polys[1] ≈ (x[1] - x[3]) / √2 #src +polys = gram[3].basis.bases[].elements #src +@test length(polys) == 1 #src +@test polys[] ≈ (x[1] - x[2] + x[3] - x[4]) / 2 #src +@test size(gram[3].Q) == (1, 1) #src +@test gram[3].Q[1, 1] ≈ 1.0 atol=1e-6 #src +gram_matrix(con_ref) diff --git a/src/Certificate/Symmetry/Symmetry.jl b/src/Certificate/Symmetry/Symmetry.jl index 9a49c76d1..c6119c288 100644 --- a/src/Certificate/Symmetry/Symmetry.jl +++ b/src/Certificate/Symmetry/Symmetry.jl @@ -3,6 +3,7 @@ module Symmetry import LinearAlgebra import MutableArithmetics as MA +import StarAlgebras as SA import MultivariatePolynomials as MP import MultivariateBases as MB diff --git a/src/Certificate/Symmetry/wedderburn.jl b/src/Certificate/Symmetry/wedderburn.jl index 059088dd8..bb50e8b66 100644 --- a/src/Certificate/Symmetry/wedderburn.jl +++ b/src/Certificate/Symmetry/wedderburn.jl @@ -1,3 +1,18 @@ +function SymbolicWedderburn.decompose( + p::MB.Polynomial, + hom::SymbolicWedderburn.InducedActionHomomorphism, +) + return [hom[p]], [1] +end + +function SymbolicWedderburn.decompose( + p::SA.AlgebraElement, + hom::SymbolicWedderburn.InducedActionHomomorphism, +) + return [hom[k] for k in SA.supp(p)], + [v for (_, v) in SA.nonzero_pairs(SA.coeffs(p))] +end + function SymbolicWedderburn.decompose( k::MP.AbstractPolynomialLike, hom::SymbolicWedderburn.InducedActionHomomorphism, @@ -10,14 +25,6 @@ function SymbolicWedderburn.decompose( return indcs, coeffs end -function SymbolicWedderburn.ExtensionHomomorphism( - action::SymbolicWedderburn.Action, - basis::MB.SubBasis{MB.Monomial}, -) - monos = collect(basis.monomials) - return SymbolicWedderburn.ExtensionHomomorphism(Int, action, monos) -end - struct VariablePermutation <: SymbolicWedderburn.ByPermutations end _map_idx(f, v::AbstractVector) = map(f, eachindex(v)) _tuple_map_idx(f, ::Tuple{}, i) = tuple() @@ -35,6 +42,18 @@ function SymbolicWedderburn.action( end abstract type OnMonomials <: SymbolicWedderburn.ByLinearTransformation end +function SymbolicWedderburn.action( + a::Union{VariablePermutation,OnMonomials}, + el, + p::MB.Polynomial{MB.Monomial}, +) + res = SymbolicWedderburn.action(a, el, p.monomial) + if res isa MP.AbstractMonomial + return MB.Polynomial{MB.Monomial}(res) + else + return MB.algebra_element(res) + end +end function SymbolicWedderburn.action( a::Union{VariablePermutation,OnMonomials}, el, @@ -75,17 +94,37 @@ end function SumOfSquares.matrix_cone_type(::Type{<:Ideal{C}}) where {C} return SumOfSquares.matrix_cone_type(C) end + +function _multi_basis_type(::Type{<:MB.SubBasis{B,M}}, ::Type{T}) where {B,M,T} + SC = SA.SparseCoefficients{M,T,Vector{M},Vector{T}} + AE = SA.AlgebraElement{MB.Algebra{MB.FullBasis{B,M},B,M},T,SC} + return Vector{MB.SemisimpleBasis{AE,Int,MB.FixedBasis{B,M,T,SC}}} +end function MA.promote_operation( ::typeof(SumOfSquares.Certificate.gram_basis), - ::Type{<:Ideal}, -) - return Vector{Vector{MB.FixedPolynomialBasis}} + C::Type{<:Ideal{SubC}}, +) where {SubC} + return _multi_basis_type( + MA.promote_operation(SumOfSquares.Certificate.gram_basis, SubC), + _coeff_type(C), + ) end function MA.promote_operation( ::typeof(SumOfSquares.Certificate.zero_basis), - ::Type{<:Ideal}, -) - return MB.Monomial + ::Type{<:Ideal{C}}, + ::Type{B}, + ::Type{D}, + ::Type{G}, + ::Type{W}, +) where {C,B,D,G,W} + return MA.promote_operation( + SumOfSquares.Certificate.zero_basis, + C, + B, + D, + G, + W, + ) end SumOfSquares.Certificate.zero_basis(::Ideal) = MB.Monomial function SumOfSquares.Certificate.reduced_polynomial( @@ -146,13 +185,19 @@ function matrix_reps(pattern, R, basis, ::Type{T}, form) where {T} end end +function _coeff_type(C::Type{<:Ideal}) + return SumOfSquares._complex(Float64, SumOfSquares.matrix_cone_type(C)) +end + function SumOfSquares.Certificate.gram_basis(cert::Ideal, poly) basis = SumOfSquares.Certificate.gram_basis(cert.certificate, poly) - T = SumOfSquares._complex( - Float64, - SumOfSquares.matrix_cone_type(typeof(cert)), - ) - return _gram_basis(cert.pattern, basis, T) + return _gram_basis(cert.pattern, basis, _coeff_type(typeof(cert))) +end + +function _fixed_basis(F, basis) + return MB.FixedBasis([ + MB.implicit(MB.algebra_element(row, basis)) for row in eachrow(F) + ]) end function _gram_basis(pattern::Pattern, basis, ::Type{T}) where {T} @@ -231,14 +276,18 @@ function _gram_basis(pattern::Pattern, basis, ::Type{T}) where {T} # Moreover, `Q = kron(C, I)` is not block diagonal but we can get a block-diagonal # `Q = kron(I, Q)` by permuting the rows and columns: # `(U[1:d:(1+d*(m-1))] * F)' * basis.monomials`, `(U[2:d:(2+d*(m-1))] * F)' * basis.monomials`, ... - map(1:d) do i - return MB.FixedPolynomialBasis( - (transpose(U[:, i:d:(i+d*(m-1))]) * F) * basis.monomials, - ) - end + return MB.SemisimpleBasis( + map(1:d) do i + return _fixed_basis( + transpose(U[:, i:d:(i+d*(m-1))]) * F, + basis, + ) + end, + ) else - F = convert(Matrix{T}, R) - [MB.FixedPolynomialBasis(F * basis.monomials)] + return MB.SemisimpleBasis([ + _fixed_basis(convert(Matrix{T}, R), basis), + ]) end end end diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index 4dcb8c32f..20bc096d3 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -16,7 +16,7 @@ Base.:+(::_NonZero, a::_NonZero) = a function _combine_with_gram( basis::MB.SubBasis{B,M}, - gram_bases::AbstractVector{<:MB.SubBasis}, + gram_bases::AbstractVector{<:SA.ExplicitBasis}, weights, ) where {B,M} p = zero(_NonZero, MB.algebra(MB.FullBasis{B,M}())) diff --git a/src/gram_matrix.jl b/src/gram_matrix.jl index 9648bce1c..64566c523 100644 --- a/src/gram_matrix.jl +++ b/src/gram_matrix.jl @@ -137,26 +137,71 @@ function GramMatrix(Q::AbstractMatrix, monos::AbstractVector) return GramMatrix(Q, MB.SubBasis{MB.Monomial}(sorted_monos), σ) end -function _term_element(α, p::MB.Polynomial{B,M}) where {B,M} - return MB.algebra_element( - MB.sparse_coefficients(MP.term(α, p.monomial)), - MB.FullBasis{B,M}(), +# `_gram_operate!` is a temporary workaround waiting for a cleaner solution +# in https://github.com/JuliaAlgebra/StarAlgebras.jl/pull/62 + +function _gram_operate!( + op, + p, + α, + row::MB.Polynomial, + col::MB.Polynomial, + args::Vararg{Any,N}, +) where {N} + return MA.operate!( + op, + p, + MB.term_element(true, row), # TODO `*` shouldn't be defined for group elements so this is a hack + MB.term_element(α, col), + args..., + ) +end + +function _gram_operate!( + op, + p, + α, + row::SA.AlgebraElement, + col::SA.AlgebraElement, + args::Vararg{Any,N}, +) where {N} + return MA.operate!( + op, + p, + true * row, # TODO `*` shouldn't be defined for group elements so this is a hack + α * col, + args..., ) end +function _gram_operate!( + op, + p, + α, + row::MB.SemisimpleElement, + col::MB.SemisimpleElement, + args::Vararg{Any,N}, +) where {N} + for (r, c) in zip(row.elements, col.elements) + _gram_operate!(op, p, α, r, c, args...) + end +end + function MA.operate!( op::SA.UnsafeAddMul{typeof(*)}, p::SA.AlgebraElement, g::GramMatrix, args::Vararg{Any,N}, ) where {N} - for col in eachindex(g.basis) - for row in eachindex(g.basis) - MA.operate!( + for row in eachindex(g.basis) + row_star = SA.star(g.basis[row]) + for col in eachindex(g.basis) + _gram_operate!( op, p, - _term_element(true, SA.star(g.basis[row])), - _term_element(g.Q[row, col], g.basis[col]), + g.Q[row, col], + row_star, + g.basis[col], args..., ) end