Simplify and fix similar and zero for our matrix types

[fingolfin]

Remove the internal helper _similar which untangles the MatElem and MatRingElem code further and makes it easier to see what's going on.

Also change similar for MatElem to default to calling zero_matrix. With this similar(::ZZMatrix) in Nemo returns a new ZZMatrix instead of a AbstractAlgebra.Generic.MatSpaceElem{ZZRingElem}.

Resolves Nemocas/Nemo.jl#1826

[lgoettgens]

From a first glance it looks like that some tests depend on the old (wrong) behaviour

[thofma]

There is also a doctests, which illustrates the "feature", that sometimes entries can be #undef.

[fingolfin]

Yeah so one thing that is wrong is that the similar(x) method delegates to similar(x::MatElem, R::NCRing, r::Int, c::Int), but that now is implemented thus:

similar(x::MatElem, R::NCRing, r::Int, c::Int) = zero_matrix(R, r, c)

So it ignores it first argument.

This is apparently even desired (?) but only if the user supplied R ?!? More on that below, but the key here is: if only a matrix is provided as input, then the resulting output matrix should have the same type. Which the current similar method didn't really ensure either in general, it just always made a Generic.MatSpaceElem... (and to fix this in general we need the construction via matType(R, r, c)).

It is unclear to me whether similar(mat, r, c) should preserve the matrix type or not.

---

So here is for some background:

Those failing tests in Nemo reference Nemocas/Nemo.jl#651 which describes:

It was decided recently (Nemocas/Nemo.jl#638) that similar(mat, ring) doesn't necessarily return a matrix of the the same type as mat.

There is the also a linked issue #490 which was also closed long ago. It contains a somewhat relevant example perhaps motivating this (slightly edited for context):

[Tommy] wants change_base_ring(ZZ, M) in Nemo, where M is a Generic.MatSpaceElem{nf_elem} to return an fmpz_mat. Reason: both fmpz_mat and Generic.MatSpaceElem{nf_elem} are matrices that Nemo would create.

[Bill wants] change_base_ring(Nemo.ZZ, p) in Singular, where p is a Singular spoly{fmpq} to return an spoly{fmpz}. Reason: both polynomials are what Singular would create over a Nemo coefficient ring.

So we don't disagree on what we want, it's just not implemented that way at present (and probably wasn't implemented correctly before the new similar functionality, either).

Of course that's for change_base_ring but my understanding is that similar(mat, R, r, c) is meant to work along the same lines (verbatim for the first example, and perhaps adjusted to smatrix for the second example?).

Though it is not quite clear to me why similar should do this... but that ship has sailed, we are stuck with it for now.

That leaves the task of how to specify / document / implement this coherently, but that goes a bit beyond this PR. I'll see what I can do here at least in the coming days/weeks.

[fingolfin]

If I understand it right then similiar(x::Generic.MatSpaceElem{ZZRingElem}) should return another Generic.MatSpaceElem{ZZRingElem} but similiar(x::Generic.MatSpaceElem{ZZRingElem}), ZZ) should return a ZZMatrix.

That then means we must not document it as similar(x::MatRingElem, R::NCRing=base_ring(x)) as that's just wrong now.

Presumably the same for similar(x::MatElem, r::Int, c::Int) (preserve type) versus similar(x::MatElem, R::NCRing, r::Int, c::Int) (may change matrix type).

For the cases taking a ring, to what matrix types shall it change? What we can do for now is to produce something of dense_matrix_type(R), that would also be easy to implement, and fits with the first example (the one involving fmpz_mat). But it would not fit with the second example; presumably there, similar(x::smatrix{ZZRingElem}, ::QQField) should return a smatrix{QQFieldElem} but that's not something I can easily explain in the AA documentation... Huhrm

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 88.16%. Comparing base (174cdc9) to head (fd156e4).
Report is 9 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1890      +/-   ##
==========================================
+ Coverage   88.13%   88.16%   +0.03%     
==========================================
  Files         120      120              
  Lines       30294    30295       +1     
==========================================
+ Hits        26700    26710      +10     
+ Misses       3594     3585       -9     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[fingolfin]

This should be ready now :-)

[thofma]

Fixes the Nemo issue, so it is an improvement.

[fingolfin]

Hrm, I just realized another twist to the rules for similar: for similar on matrix views, we want the result to not be a view anymore (I think, based on what we do right now and what Julia does, and because it seems the one sensibel thing to do)

But that means the type of the resulting matrix differs from the input matrix even in the single-argument version similar(myMatrixView).

Can be dealt with of course, but it's yet another wrinkle in the story that our docs don't tell (I think? Didn't search thoroughly yet, so perhaps it is better than I claim)

[fingolfin]

A very simple rule would be if we ignore the smatrix example by Bill, and just use these rules (again):

• similar(mat) = similar(mat, base_ring(mat))
• the result of similar(mat, R) is always of type dense_matrix_type(R)

But then there are trouble again with Generic.MatSpaceElem{ZZRingElem} e.g. then -mat suddenly has a different type then mat.

Urgh

[thofma]

Maybe a "proper" fix is to have AbstractAlgebra, Nemo and Singular have their own similar. Another fix would be to remove similar altogether.

[fingolfin]

Also, let's not throw out the baby with the wash water :-).

First off, I think this PR here still is "safe" (and in any case, even if something regressed, it should be easy enough to adjust). So I'll merge it.

Next, I note that also similar in Julia is not 100% clear, but overall I find its definition workable:

  similar(array, [element_type=eltype(array)], [dims=size(array)])

  Create an uninitialized mutable array with the given element type and size, based upon
  the given source array. The second and third arguments are both optional, defaulting to
  the given array's eltype and size. The dimensions may be specified either as a single
  tuple argument or as a series of integer arguments.

  Custom AbstractArray subtypes may choose which specific array type is best-suited to
  return for the given element type and dimensionality. If they do not specialize this
  method, the default is an Array{element_type}(undef, dims...).

What really causes me confusion is the different between similar(mat) and similar(mat, base_ring(R)) -- that they can return different values. Indeed whoever wrote the documentation also did not consider this as they were documented to return the same thing.

To me it feels like this is what we need to resolve. Perhaps by indeed splitting similar into (at least) two functions.

1. one function which maps a Generic.MatSpaceElem{ZZRingElem}) to another -- i.e. preserving the type here
   • this function would still be allowed to map a view into something else
   • this is presumably mainly useful for internal constructions e.g. when implementing arithmetic like -mat ?
2. one function which maps a Generic.MatSpaceElem{ZZRingElem}) to a ZZMatrix -- i.e. choosing an "optimal" type
   • this is probably what we want "most of the time" (?)
   • and the type we'd use would be dense_matrix_type(R) which is also easy to explain and document.

In each case, similar(mat) and similar(mat, base_ring(R)) would do the same thing.

My feeling is that we want 2, and we may not even need 1 -- adjusting the artihmetic code in e.g. Generic.MatSpaceElem to not use similar is easy enough. For other matrix types, well, they "know" whether they are "the" blessed dense_matrix_type and then of course can use similar if they want to.

[fingolfin]

Maybe we should just discuss the desired behavior in triage?

[thofma]

I don't think a discussion of this during the triage meeting would be productive.

For the discussion, I think we should not assume that something like Generic.MatSpaceElem{ZZRingElem} can exist, but that all matrices over R have type dense_matrix_type(R). I think we then should always have

typeof(similar(M::dense_matrix_type(R))) == dense_matrix_type(R)

and

typeof(similar(M, R)) == dense_matrix_type(R)

The philosophy is that for a user there is only one matrix type per element type and this is preserved under similar.

If AA or Singular really need a method that makes similar on smatrix or Generic.MatSpaceElem preserve the type, we can add it.

[lgoettgens]

I agree with @thofma.
But I think adding a more general dispatch for all MatElem that satisfies
typeof(similar(M)) = typeof(M)
Would be beneficial for internal arithmetics of AA. (eg for doing similar inside the definition of -(::Generic.MatElem)

[fingolfin]

I think we agree so far that

typeof(similar(M::dense_matrix_type(R))) == dense_matrix_type(R)
typeof(similar(M, R)) == dense_matrix_type(R)

should hold. Good!

That leaves the question whether to require

1. typeof(similar(M)) = typeof(M) holds unless M is a view, or
2. typeof(similar(M)) = dense_matrix_type(base_ring(M)) holds

But this "choice" only matters because we allow instances of Generic.MatSpaceElem{ZZRingElem} to exist. I think @thofma hit the core of the issue there -- if we simply don't allow such instances to exist, the problem mostly goes away, as then both 1 and 2 can be true simultaneously (well, sort off -- I am ignoring "special" cases like smatrix here, but I think that's fine. Also, I would document it with rule 2 as that is far simpler and does not require understanding what "view" is etc.).

Put concretely, we could change the Generic.MatSpaceElem constructor to throw an error if one tries to instantiate a Generic.MatSpaceElem{T} for which dense_matrix_type(T) != Generic.MatSpaceElem{T}.

Of course this would be mildly breaking (as in: some code might accidentally have some of those; and the Nemo test suite on purpose tries to create such object). But other than that I can't think of any serious downside?

[fingolfin]

This method would go if we agree on the proposal to forbid MatSpaceElem{T} instances if dense_matrix_type(T) != MatSpaceElem{T}.

But in the meantime it is necessary to preserve Nemo tests.

[fingolfin]

To clarify: this PR does not implement the similar revision we are discussing. I've now adjust it once more to cleanup some docstrings. The docstrings are currently not quite correct. But I deliberate leave them as they are, as this PR has dragged on for far too long already.

If @lgoettgens can get himself to agree, I suggest we merge this now as it improves the state of affairs (I've suggested that before but then didn't merge... ah well).

For the docstrings, I note that also the manual itself currently says many no longer quite true things about matrices and matrix ring elements. I would adjust these and the docstrings in one go. But it doesn't make sense to adjust them to the current state of things if then later agree to implement changes to how similar works. Hence my reluctance for further changes in this PR.

If we then agree on the plan to "forbid Generic.MatSpaceElem{ZZRingElem}" then in fact the docstrings here are mostly fine, I would just extend them to spell out more clearly what is going on.

[fieker]

On Tue, Nov 12, 2024 at 12:39:50AM -0800, Max Horn wrote: To clarify: this PR does not implement the `similar` revision we are discussing. I've now adjust it once more to cleanup some docstrings. The docstrings are currently not quite correct. But I deliberate leave them as they are, as this PR has dragged on for far too long already. If @lgoettgens can get himself to agree, I suggest we merge this now as it improves the state of affairs (I've suggested that before but then didn't merge... ah well). For the docstrings, I note that also the manual itself currently says many no longer quite true things about matrices and matrix ring elements. I would adjust these and the docstrings in one go. But it doesn't make sense to adjust them to the current state of things if then later agree to implement changes to how `similar` works. Hence my reluctance for further changes *in this PR*. If we then agree on the plan to "forbid `Generic.MatSpaceElem{ZZRingElem}`" then in fact the docstrings here are mostly fine, I would just extend them to spell out more clearly what is going on.

I don't thinkn we're going to put effort in to forbidding it - but we don't care about it and will not support it (actively)

…

-- Reply to this email directly or view it on GitHub: #1890 (comment) You are receiving this because you are subscribed to this thread. Message ID: ***@***.***>

[thofma]

I agree with your summary @fingolfin, but also agree with @fieker that we probably should not enforce this by throwing an error. I can imagine that for some horrible internal hacks it might be advantage to still have this option. (Similar to how we allow Generic.PolyRing{ZZRingElem}).

[lgoettgens]

This could be changed in the future to use the UndefInitializer, right? (but not needed right now)
so something like dense_matrix_type(R)(R, undef, r, c)

[fingolfin]

Yes, that's the plan (but we can't change it now because that'd be a breaking change and wouldn't work with older Nemo versions.