Nemocas / AbstractAlgebra.jl

Generic abstract algebra functionality in pure Julia (no C dependencies)
https://nemocas.github.io/AbstractAlgebra.jl/dev/index.html
Other
161 stars 62 forks source link

Can we multiply a vector by a matrix? #572

Open ulthiel opened 4 years ago

ulthiel commented 4 years ago

The following seems natural to me but doesn't work:

V=VectorSpace(QQ,3)
A=matrix(QQ,3,3,[0,1,-1, 1,0,-1, 0,0,-1])
gen(V,1)*A

Let's be able to actually act on V!

fieker commented 4 years ago

On Thu, Feb 06, 2020 at 04:49:07AM -0800, Ulrich Thiel wrote:

The following seems natural to me but doesn't work:

V=VectorSpace(QQ,3)
A=matrix(QQ,3,3,[0,1,-1, 1,0,-1, 0,0,-1])
gen(V,1)*A

Let's be able to actually act on V!

On thinking about this: cannot made to work What you want is to turn A into an automorphism (a map) from V -> V otherwise the product will be a matrix. There is no reason why the product is in V. ... unless A is actually an automorphism...

-- You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/572

ulthiel commented 4 years ago

Why does A need to be an automorphism? It defines a map V -> V (V has a fixed basis).

In Magma

> V:=VectorSpace(Rationals(),3);
> A:=Matrix(Rationals(),3,3,[0,1,-1, 1,0,-1, 0,0,-1]);
> v := V.1;
> v*A;
( 0  1 -1)
thofma commented 4 years ago

Since VectorSpace(F, n) always produces a free modulo of rank n, we have a fixed basis and there is no ambiguity.

We can just make it work for any ring.

ulthiel commented 4 years ago

Since VectorSpace(F, n) always produces a free modulo of rank n, we have a fixed basis and there is no ambiguity.

We can just make it work for any ring.

Is this supporting my comment or Claus's?

thofma commented 4 years ago

Is this supporting my comment or Claus's?

Sorry, our posts crossed. This was meant to support your proposal. I don't think there is a good reason to not support it.

fieker commented 4 years ago

On Thu, Feb 06, 2020 at 05:50:10AM -0800, thofma wrote:

Is this supporting my comment or Claus's?

Sorry, our posts crossed. This was meant to support your proposal. I don't think there is a good reason to not support it.

julia> Q = FreeModule(QQ, 3, cached = false) julia> R = FreeModule(QQ, 3, cached = false) julia> m = matrix(...) gen(Q[1])*m in Q? or in R?

The rule being that any n x n matrix acts as on any n-dim vectorspace?

Here Q !== R

-- You are receiving this because you commented. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/572#issuecomment-582913945

thofma commented 4 years ago

The result will be in Q. There is a canonical action of Mat(F, n, n) on every F^n. I don't think this is a problem.

It is the same as for the other canonical promotions/embeddings that we automatically do: The polynomial 2 * x will be an element of Qx and not Qy, where

Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
Qy, y = PolynomialRing(FlintQQ, "y", cached = false)
fieker commented 4 years ago

On Thu, Feb 06, 2020 at 06:01:53AM -0800, thofma wrote:

The result will be in Q. There is a canonical action of Mat(F, n, n) on every F^n. I don't think this is a problem.

It is the same as for the other canonical promotions/embeddings that we automatically do: The polynomial 2 * x will be an element of Qx and not Qy, where

Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
Qy, y = PolynomialRing(FlintQQ, "y", cached = false)

OK, so square matrices (only square ones!) act as automorphisms on any VS of the correct dim (over the same field)? -- You are receiving this because you commented. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/572#issuecomment-582918633

ulthiel commented 4 years ago

So, I'm supposed to do

f = ModuleHomomorphism(V, V, A)
f(v)

Okay, makes sense in a way. But I feel that many users will fall into this trap; my first attempt is how you do it in basically every other system I would say. If you want it like this then maybe a "common pitfalls" section in the documentation would be helpful?

thofma commented 4 years ago

Yes, we should probably document this for the time being.

Now I know what you mean @fieker. You should have mentioned the non-square case right at the beginning. Does Magma have a unique object VectorSpace(F, n) for every pair (F, n)?

ulthiel commented 4 years ago

Sorry to interrupt but since we're at it: how do I get back the matrix of a vector space homomorphism as above? matrix(f) doesn't work. Also, inv(f) doesn't work. Do I really have to create it as ModuleIsomorphism before I can do inv(f)?

thofma commented 4 years ago

@wbhart Is there a clean way to do this with the current interface?

wbhart commented 4 years ago

Computationally, I think there is a difference between:

We don't distinguish these things at present. Currently there is only one free vector space of a given dimension over a given field in the system, unless you turn caching off

This implies that currently it would be safe to assume that a matrix is a linear map between the unique vector spaces over its coefficient ring implied by its dimensions. In other words, I see no serious problem with @ulthiel 's proposal.... at present.

However, this violates one of the guiding principles of Oscar (as enunciated to me by others): objects should not have a mathematical meaning that changes depending on what way you view them.

We've had exactly this discussion about how in a certain well-known system for commutative algebra, a matrix is a homomorphism, and how we aren't going to make this same "mistake" in Oscar.

So if we go down this road, we need to understand the mathematical consequences. The most obvious one is that we encourage users to implement algorithms with a set of computational tools (like matrices), instead of mathematical tools (like homomorphisms).

One possibility is to separate these two worlds. The highest level of Oscar could present things in well-defined mathematical concepts. The libraries (Nemo/Singular/AbstractAlgebra, etc) could present computational tools to build Oscar on top of.

Of course the risk is no one wants to use Oscar because it is too particular and no one wants to use the libraries because they aren't mathematical enough.

I do feel that a lot of serious mathematicians don't get involved in computer algebra because they feel the systems aren't mathematical enough, or alternatively because they try to use a computer algebra system and discover it doesn't abuse notation, make identifications and use the conventions that they are used to.

In answer to @thofma , yeah it can be done cleanly enough. You have to construct the parent of the output, but that's not computationally demanding. It's just a dictionary lookup based on ring and dimension.

I made the suggestion that we optionally only support this when caching is on. But I don't see a clean way to do this, other than demand that someone already constructed the output ring (so that it's in the cache). I don't like the idea of explaining to users why this is though.

One thing we can certainly agree on is that if we don't support it, we should make it absolutely trivial to construct the homomorphism from the matrix. And perhaps we could give a more meaningful error message, e.g. "higher dimensional tensors not implemented", or "field does not possess the ABDF property", or something equally mystifying so that the user knows that whatever they thought this should do, someone else thought it meant something else. (Just joking of course.)

wbhart commented 4 years ago

In summary, I'm certainly ok with supporting @ulthiel 's proposal in AbstractAlgebra/Nemo. However, I don't see how it doesn't leak into Oscar. Operators, like * are handled differently than functions. As far as I know, we can't import only the definitions we want into Oscar.

Therefore (IMHO), this has to go to a meeting for a decision, with all parties present who have a stake in the decision.

thofma commented 4 years ago

We should just implement it here in AA, even for the non-square case. People working with things where the caching is turned off cannot use the syntax v * A anyway. So not supporting it does not help anyone, but just inconveniences the ordinary user.

P.S.: I don't think that operators like * are handled differently than functions. I think they are literally just functions.

fingolfin commented 3 years ago

Just to say, this functionality is already in Oscar.jl (and has been for some time), in src/Groups/matrices/matrix_manipulation.jl; this slipped past me in review, I only recently noticed it. Here's the code -- I'll work on removing at least some of them (e.g. the getindex), but the ability to apply a matrix (or even a matrix group element) to a vector seems like something that a lot of people would want to do, and for which I don't see the harm it'd cause...?

Base.getindex(V::AbstractAlgebra.Generic.FreeModule, i::Int) = gen(V, i)

# scalar product
Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = (v.v*transpose(u.v))[1]

Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},x::MatElem{T}) where T <: FieldElem = v.parent(v.v*x)
Base.:*(x::MatElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = x*transpose(u.v)

# evaluation of the form x into the vectors v and u
Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},x::MatElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = (v.v*x*transpose(u.v))[1]

Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},x::MatrixGroupElem{T}) where T <: FieldElem = v.parent(v.v*x.elm)
Base.:*(x::MatrixGroupElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = x.elm*transpose(u.v)

# evaluation of the form x into the vectors v and u
Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},x::MatrixGroupElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = (v.v*x.elm*transpose(u.v))[1]

map(f::Function, v::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = v.parent(map(f,v.v))
fieker commented 3 years ago

In the end, I don't mind too much either way. Problem that I see are e.g. kernel(mat) => mat kernel(morphism) => subspace same for image What do we want here? I guess at the core is the ambiguity of vectors:

I see the convenience - and the range of future and current problems.

On Thu, Aug 26, 2021 at 05:47:45PM -0700, Max Horn wrote:

Just to say, this functionality is already in Oscar.jl (and has been for some time), in src/Groups/matrices/matrix_manipulation.jl; this slipped past me in review, I only recently noticed it. Here's the code -- I'll work on removing at least some of them (e.g. the getindex), but the ability to apply a matrix (or even a matrix group element) to a vector seems like something that a lot of people would want to do, and for which I don't see the harm it'd cause...?

Base.getindex(V::AbstractAlgebra.Generic.FreeModule, i::Int) = gen(V, i)

# scalar product
Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = (v.v*transpose(u.v))[1]

Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},x::MatElem{T}) where T <: FieldElem = v.parent(v.v*x)
Base.:*(x::MatElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = x*transpose(u.v)

# evaluation of the form x into the vectors v and u
Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},x::MatElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = (v.v*x*transpose(u.v))[1]

Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},x::MatrixGroupElem{T}) where T <: FieldElem = v.parent(v.v*x.elm)
Base.:*(x::MatrixGroupElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = x.elm*transpose(u.v)

# evaluation of the form x into the vectors v and u
Base.:*(v::AbstractAlgebra.Generic.FreeModuleElem{T},x::MatrixGroupElem{T},u::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = (v.v*x.elm*transpose(u.v))[1]

map(f::Function, v::AbstractAlgebra.Generic.FreeModuleElem{T}) where T <: FieldElem = v.parent(map(f,v.v))

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/572#issuecomment-906838660

thofma commented 3 years ago

As I wrote above, I am happy to add this. (Some of the functions added in Oscar are wrong for non-square matrices and the v * A * u should also be removed.)

Since most people will not work with FreeModuleElem but with Vector, the question is what to do about

    *(v::Vector, M::MatrixElem)

so, Base Vector times Oscar matrix. As expected, for Base a Vector is a column vector and is hence multiplied with a Base Matrix from the right. Of course this clashes with our notion of a vector, since vectors are for us rows. Should we allow v::Vector * M::MatrixElem or M::MatrixElem * v::Vector?

Or we just allow v * M and M * v by letting matrices act on F^n on both sides. Does this lead to some inconsistencies?

fieker commented 3 years ago

On Fri, Aug 27, 2021 at 01:32:16AM -0700, Tommy Hofmann wrote: I'd allow Mat * Vec as this is the math definition. I am sorry that julia thinks in columns...

As I wrote above, I am happy to add this. (Some of the functions added in Oscar are wrong for non-square matrices and the v * A * u should also be removed.)

Since most people will not work with FreeModuleElem but with Vector, the question is what to do about

    *(v::Vector, M::MatrixElem)

so, Base Vector times Oscar matrix. As expected, for Base a Vector is a column vector and is hence multiplied with a Base Matrix from the right. Of course this clashes with our notion of a vector, since vectors are for us rows. Should we allow v::Vector * M::MatrixElem or M::MatrixElem * v::Vector?

Or we just allow v * M and M * v by letting matrices act on F^n on both sides. Does this lead to some inconsistencies?

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/572#issuecomment-907028069

tthsqe12 commented 3 years ago

2 cents: when I was putting these in flint, I had the idea of using Julia's Vectors, since we don't have to think about what the correct parent is for these Vectors. If you have a matrix with entry type T and a Vector{T}, you should be able to multiply them (in either order) and get a Vector{T}. https://github.com/wbhart/flint2/pull/961