Rework linear independence

src/DeformationBases/ArcDiagDeformBasis.jl

@@ -102,8 +102,8 @@ struct ArcDiagDeformBasis{C <: RingElem} <: DeformBasis{C}
         if !no_normalize
             iter = unique(Iterators.filter(b -> !iszero(b), iter))
             collected = Vector{DeformationMap{C}}(collect(iter))
-            _, rels = is_linearly_independent_with_relations(coefficient_ring(sp), reverse(collected))
-            inds = [1 + ncols(rels) - (findfirst(!iszero, vec(rels[i, :]))::Int) for i in nrows(rels):-1:1]
+            _, rels = is_linearly_independent_with_relations(coefficient_ring(sp), collected)
+            inds = [findlast(!iszero, vec(rels[i, :]))::Int for i in 1:nrows(rels)]
             deleteat!(collected, inds)
             return new{C}(length(collected), collected, extra_data, normalize)
         end

src/LinearIndependence.jl

@@ -1,71 +1,100 @@
-function is_linearly_independent(F::Field, V::Vector{T}) where {T}
-    return is_linearly_independent_with_relations(F, V)[1]
+const _linear_independence_rref_cutoff = 0
+
+function column_rref!(mat::MatElem{T}) where {T <: FieldElem}
+    rk = rref!(AbstractAlgebra.Solve.lazy_transpose(mat))
+    return view(mat, :, 1:rk)
+end
+
+function column_rref!(mat::Nemo._MatTypes) # change to Nemo._FieldMatTypes
+    trmat = transpose!(mat)
+    rk = rref!(trmat)
+    return view(transpose!(trmat), :, 1:rk)
 end
 
 function is_linearly_independent(V::Vector{T}) where {T}
     return is_linearly_independent_with_relations(V)[1]
 end
 
-function is_linearly_independent_with_relations(F::Field, V::Vector{<:FieldElem})
-    n, M = kernel(matrix(F, length(V), 1, V); side=:left)
-    return n == 0, M
+function is_linearly_independent(F::Field, V::Vector{T}) where {T}
+    return is_linearly_independent_with_relations(F, V)[1]
 end
 
-function is_linearly_independent_with_relations(F::Field, V::Vector{<:PolyRingElem})
-    n = length(V)
-    if n == 0
-        return false, zero_matrix(F, 0, 0)
-    end
-    rels = identity_matrix(F, n)
-    maxdeg = maximum(degree, V)
-    for i in 1:maxdeg
-        rels = basis_intersection(rels, is_linearly_independent_with_relations(F, [coeff(v, i) for v in V])[2])
-        iszero(rels) && return true, rels
-    end
-    return false, rels
+function is_linearly_independent_with_relations(V::Vector{T}) where {T}
+    M = kernel(_linear_independence_coeff_matrix(V); side=:left)
+    return nrows(M) == 0, M
+end
+
+function is_linearly_independent_with_relations(F::Field, V::Vector{T}) where {T}
+    M = kernel(_linear_independence_coeff_matrix(F, V); side=:left)
+    return nrows(M) == 0, M
 end
 
-function is_linearly_independent_with_relations(F::Field, V::Vector{<:MatElem})
+function _linear_independence_coeff_matrix(V::Vector{<:FieldElem})
+    @req length(V) > 0 "For empty vectors, the field needs to be specified"
+    return _linear_independence_coeff_matrix(parent(V[1]), V)
+end
+
+function _linear_independence_coeff_matrix(F::Field, V::Vector{<:FieldElem})
+    return matrix(F, length(V), 1, V)
+end
+
+function _linear_independence_coeff_matrix(F::Field, V::Vector{<:PolyRingElem})
+    @req isempty(V) || all(v -> parent(v) === parent(V[1]), V) "Incompatible polynomial rings"
+    return reduce(
+        hcat,
+        begin
+            mat = _linear_independence_coeff_matrix(F, [coeff(v, i) for v in V])
+            if ncols(mat) > _linear_independence_rref_cutoff * nrows(mat)
+                column_rref!(mat)
+            else
+                mat
+            end
+        end for i in 0:maximum(degree, V; init=-1);
+        init=zero_matrix(F, length(V), 0),
+    )
+end
+
+function _linear_independence_coeff_matrix(F::Field, V::Vector{<:MatElem})
     n = length(V)
     if n == 0
-        return false, zero_matrix(F, 0, 0)
-    end
-    @req all(size(v) == size(V[1]) for v in V) "Size mismatch"
-    rels = identity_matrix(F, n)
-    for i in eachindex(V[1])
-        rels = basis_intersection(rels, is_linearly_independent_with_relations(F, [v[i] for v in V])[2])
-        iszero(rels) && return true, rels
+        return zero_matrix(F, n, 0)
     end
-    return false, rels
+    @req all(v -> parent(v) === parent(V[1]), V) "Incompatible matrix spaces"
+    return reduce(
+        hcat,
+        begin
+            mat = _linear_independence_coeff_matrix(F, [v[i] for v in V])
+            if ncols(mat) > _linear_independence_rref_cutoff * nrows(mat)
+                column_rref!(mat)
+            else
+                mat
+            end
+        end for i in eachindex(V[1]);
+        init=zero_matrix(F, n, 0),
+    )
 end
 
-function is_linearly_independent_with_relations(F::Field, V::Vector{<:FreeAssAlgElem})
-    coeff_maps = [Dict{Vector{Int}, elem_type(F)}(zip(exponent_words(v), coefficients(v))) for v in V]
-    words = reduce(union, keys.(coeff_maps))
+function _linear_independence_coeff_matrix(F::Field, V::Vector{<:FreeAssAlgElem})
     n = length(V)
-    rels = identity_matrix(F, n)
-    for word in words
-        rels = basis_intersection(
-            rels,
-            is_linearly_independent_with_relations(F, [get(coeff_map, word, zero(F)) for coeff_map in coeff_maps])[2],
-        )
-        iszero(rels) && return true, rels
+    if n == 0
+        return zero_matrix(F, n, 0)
     end
-    return false, rels
+    R = parent(V[1])
+    C = base_ring(R)
+    @req all(v -> parent(v) === R, V) "Incompatible algebras"
+    coeff_maps = [Dict{Vector{Int}, elem_type(F)}(zip(exponent_words(v), coefficients(v))) for v in V]
+    support_words = reduce(union, keys.(coeff_maps))
+    return reduce(
+        hcat,
+        begin
+            mat = _linear_independence_coeff_matrix(F, [get(coeff_map, word, zero(C)) for coeff_map in coeff_maps])
+            if ncols(mat) > _linear_independence_rref_cutoff * nrows(mat)
+                column_rref!(mat)
+            else
+                mat
+            end
+        end for word in support_words;
+        init=zero_matrix(F, n, 0),
+    )
 end
 
-function basis_intersection(V::MatElem{T}, W::MatElem{T}) where {T <: FieldElem}
-    # using Zassenhaus
-    @req base_ring(V) == base_ring(W) "Incompatible base rings"
-    @req ncols(V) == ncols(W) "Size mismatch"
-    if nrows(V) == 0 || nrows(W) == 0
-        return zero(V, 0, ncols(V))
-    end
-    if V == W
-        return V
-    end
-    block_mat = vcat(hcat(V, V), hcat(W, zero(W)))
-    s = rref!(block_mat)
-    r = rank(block_mat[:, 1:ncols(V)])
-    return block_mat[r+1:s, ncols(V)+1:end]
-end