Nemocas / Nemo.jl

Julia bindings for various mathematical libraries (including flint2)
http://nemocas.github.io/Nemo.jl/
Other
176 stars 57 forks source link

Matrix times transpose error #1816

Closed linus-md closed 1 day ago

linus-md commented 2 days ago

Hello,

I want to compute A*A^T where a is a matrix of polynomial elements like that

using Nemo
n, m = 2, 2
R, A, B =  polynomial_ring(GF(101),:A=>(1:n, 1:m), :B=>(1:n, 1:m))
A * transpose(A)
ERROR: MethodError: no method matching copy(::FqMPolyRingElem)

which does not work, why is that?

1-element ExceptionStack:
MethodError: no method matching copy(::FqMPolyRingElem)

Closest candidates are:
  copy(::Random.Xoshiro)
   @ Random /usr/local/Cellar/julia/1.10.4/share/julia/stdlib/v1.10/Random/src/Xoshiro.jl:73
  copy(::Base.JuliaSyntax.SyntaxData)
   @ Base /private/tmp/julia-20240605-9481-7mnuft/julia-1.10.4/base/JuliaSyntax/src/syntax_tree.jl:217
  copy(::Markdown.MD)
   @ Markdown /usr/local/Cellar/julia/1.10.4/share/julia/stdlib/v1.10/Markdown/src/parse/parse.jl:30
  ...

Stacktrace:
 [1] matmul2x2!(C::Matrix{FqMPolyRingElem}, tA::Char, tB::Char, A::Matrix{FqMPolyRingElem}, B::Matrix{FqMPolyRingElem}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
   @ LinearAlgebra /usr/local/Cellar/julia/1.10.4/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:1003
 [2] generic_matmatmul!(C::Matrix{FqMPolyRingElem}, tA::Char, tB::Char, A::Matrix{FqMPolyRingElem}, B::Matrix{FqMPolyRingElem}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
   @ LinearAlgebra /usr/local/Cellar/julia/1.10.4/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:778
 [3] mul!
   @ /usr/local/Cellar/julia/1.10.4/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:263 [inlined]
 [4] mul!
   @ /usr/local/Cellar/julia/1.10.4/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:237 [inlined]
 [5] *(A::Matrix{FqMPolyRingElem}, B::LinearAlgebra.Transpose{FqMPolyRingElem, Matrix{FqMPolyRingElem}})
   @ LinearAlgebra /usr/local/Cellar/julia/1.10.4/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:106
 [6] top-level scope
   @ REPL[18]:1
linus-md commented 2 days ago

I guess this relates to https://github.com/Nemocas/Nemo.jl/blob/79c9770d112d106c0d2911af095b929b0b7abf6f/docs/src/developer/topics.md?plain=1#L333 this somehow but how should I mitigate the issue?

linus-md commented 1 day ago

http://nemocas.github.io/Nemo.jl/v0.6/matrix/ matrix space solve the problem!

fingolfin commented 1 day ago

You are trying to work with a Julia Matrix containing Nemo elements, something we generally don't do and hence don't test for.

The Nemo resp. Oscar way of doing this would be this:

julia> A2 = matrix(R, A)
[A[1, 1]   A[1, 2]]
[A[2, 1]   A[2, 2]]

julia> A2 * transpose(A2)
[            A[1, 1]^2 + A[1, 2]^2   A[1, 1]*A[2, 1] + A[1, 2]*A[2, 2]]
[A[1, 1]*A[2, 1] + A[1, 2]*A[2, 2]               A[2, 1]^2 + A[2, 2]^2]

That said, your code could be made to work by adding this definition:

Base.copy(x::Nemo.NCRingElem) = x

based on this method definition from the Julia stdlib:

copy(x::Number) = x # some code treats numbers as collection-like

With that added you'll get

julia> A * transpose(A)
2×2 Matrix{FqMPolyRingElem}:
 A[1, 1]^2 + A[1, 2]^2              A[1, 1]*A[2, 1] + A[1, 2]*A[2, 2]
 A[1, 1]*A[2, 1] + A[1, 2]*A[2, 2]  A[2, 1]^2 + A[2, 2]^2

I am not sure such a copy method would be safe in general, though... Seems we already have three and they either use deepcopy or create a new ring element:

src/antic/nf_elem.jl:225:Base.copy(d::AbsSimpleNumFieldElem) = deepcopy(d)
src/flint/fmpq_poly.jl:76:Base.copy(f::QQPolyRingElem) = parent(f)(f)
src/flint/fmpz.jl:201:Base.copy(a::ZZRingElem) = deepcopy(a)

Thing is, the "documentation" for Base.copy is really vague, and it does not really specify the expected semantics or what the result of it is used for :-/. In the specific case you encountered here, copy is invoked by the LinearAlgebra stdlib function matmul2x2! on matrix entries... And unfortunately it says nothing about why it does that. Looking at the LinearAlgebra documentation suggests that perhaps this is for dealing with nested matrices (matrices with matrices as entries)... Anyway, for that particular caller of copy, the method I suggested above should be fine. But there might be others were one really should be making an actual copy?

I guess we could add something like Base.copy(x::Nemo.NCRingElem) = deepcopy(x) to Nemo to be on the safe side and still allow code like yours, but the result would be even more inefficient because of all the pointless copies being made.

fingolfin commented 1 day ago

You don't need matrix spaces to solve the issue -- and be careful to not really on the documentation for the truly ancient Nemo 0.6 (from 2017) while working with current 0.45 (from 2024). The current docs can be found at http://nemocas.github.io/Nemo.jl/stable/matrix/ (for the chapter you were referencing)