`decompose` does not necessarily return a minimum-norm decomposition

[thchr]

Reading a recent preprint, I realized that the docstring for decompose, which claims it returns a minimum-norm solution, is not correct. It seems that even though the argument works for real-coefficient solutions, it does not work the solution space is restricted to integer coefficients. My previous thinking was described here: https://github.com/thchr/SymmetryBases.jl/blob/7836d9efd5ed6301ecf82c0c32e9bcf3a6a26d5d/src/symvec.jl#L374-L379

Below is an example that demonstrates that it is possible to obtain a smaller-norm solution:

using Crystalline, SymmetryBases
using LinearAlgebra
using Test

brs = bandreps(224)
n = -brs[1] + brs[6] + 2brs[9] - brs[12]
n′ = -brs[1] + brs[9] + brs[15]
@test n == n′

B = matrix(brs)
F = smith(B)
Nⁱʳʳ, Nᴱᴮᴿ = size(F.S, 1), size(F.T, 1)
Λᵍ = diagm(Nⁱʳʳ, Nᴱᴮᴿ, map(λ -> iszero(λ) ? 0.0 : inv(λ), F.SNF))'
Bᵍ = F.Tinv * Λᵍ * F.Sinv

c = Int.(Bᵍ*n) # particular solution (what we believed was the minimum-norm solution) (need not be integer-valued in general, but is here)
@test c == decompose(n, brs; allow_nonphysical=true)

# ------------------------------------------
# all solutions have the form `c + N*y` where y is an integer vector
N = round.(Int, (I-Bᵍ*B)) # null-space matrix
idxs = findall(!iszero, eachcol(N)) 
N′ = N[:,idxs] # nonzero columns of null-space

# ------------------------------------------
# Solve a mixed-integer quadratic problem (MIQP); apparently related to a broader class of
# so-called MISOCP (cf. https://discourse.julialang.org/t/constraint-seems-linear-but-the-compiler-determined-it-was-quadratic-not-sure-why-or-how-to-fix/84426/8)

using JuMP, Pajarito, HiGHS, Hypatia
m = Model( 
        # adapted from from Pajarito's README.md
        optimizer_with_attributes(
            Pajarito.Optimizer,
            "oa_solver" => optimizer_with_attributes(
                HiGHS.Optimizer,
                MOI.Silent() => true,
                "mip_feasibility_tolerance" => 1e-8,
                "mip_rel_gap" => 1e-6),
            "conic_solver" => optimizer_with_attributes(
                Hypatia.Optimizer,
                MOI.Silent() => true),
            )
)

@variable(m, v[1:length(idxs)], Int)
@objective(m, Min, sum(abs2, c+N′*v))
optimize!(m)

# obtain new coefficient solution (c′)
v_vals = round.(Int, value.(m[:v]))
c_null = N′*v_vals
c′ = c + c_null

@test B*c′ == n # true
@show norm(c)   # = 2.6457513110645907 = √7
@show norm(c′)  # = 1.7320508075688772 = √3 (new solution has smaller norm!)

I don't know what the best way to fix this robustly and yet cheaply could be. There is rarely a need to guarantee having the minimum-norm solution, but right now, the docs for decompose are not factual. Specifically, we currently, misleadingly, write: https://github.com/thchr/SymmetryBases.jl/blob/7836d9efd5ed6301ecf82c0c32e9bcf3a6a26d5d/src/symvec.jl#L340-L342

It could potentially be fixed by using the Pajarito solver above, but I don't know if the approach guarantees a solution - and it is also not super fast. Conceivably, it could be very slow in some situations.

[thchr]

Should be fixed by https://github.com/thchr/SymmetryBases.jl/commit/3f3bf09fa067a893ebcd53f85d24dc0e05da3e3d now.

[thchr]

(But beware that the solution is not unique in general: it is just some minimum norm solution - not even sure if necessarily a global minimum)