Matrices of matrices

[Waldek Hebisch]

In mail from 4 Nov 2019 Johannes Grabmeier proposed extending
two dimensional arrays to work as matrices of matrices.

I have looked at possible ways to obtain desired functionality.
One is to generalize Matrix to allow more general entries and
define arithmetic operations when needed operations are defined
in base domain. Another way is to define domain for "infinite"
matrices having only finitely many nonzero entries. Such
matrices with entries from a ring again form a ring (without unit,
that is our Rng), so "infinite" matrices would be acceptable as
base ring. For specific purpose another possibility is to
create domain of "block" matrices.

I have now mostly implemented the first possiblity, that is
generalization of Matrix. AFAICS generalized domain works
fine. However, it uncovered troubles with interpeter and
Spad compiler. Namely, Spad compiler has trouble with
condition like

if UP has zero? : UP -> Boolean then

where UP is domain parameter which satisfies
UnivariatePolynomialCategory(F) and F is UniqueFactorizationDomain.
This imples that conditon above is satisfies, but Spad
compiler can not infer this, so I had to add extra
conditionals as workarounds. There were few similar cases
(involwing different operations). Spad compiler seem
to be managable and with workarounds build goes on OK.

Interpreter problems are more troubling. First, after generalization
interpreter can generate Matrix(Symbol), which alone is positive
fact. Interpretes stronly prefers Variable or polynomials
compared to Symbol, but sometimes generates Matrix(Symbol)
in cases where previously we got Matrix(Polynomial(Integer))
or some other "mathematical" domain. Trouble is, that
when trying to do arithmetic on Matrix(Symbol) interpeter
needs coercions and coerce code is buggy. For example,
during coercion interpreter invented SQMATRIX(2, Symbol)
which is invalid type (SQMATRIX still requires base
domain satisfying Join(SemiRng, AbelianMonoid)). And
interpeter tries to use multiplcation form such illegal
domain (probably because multiplcation from SQMATRIX
is exported unconditionally) leading to "System error".
In effect we get bunch of testsuite failures. In several
other tests we get different types, like
Matrix(PositiveInteger), which is sensible but may cause
confusion. Another example of confusions is message:

The argument to the side-effect producing operation setelt! is not
allowed to be converted from type Matrix(SquareMatrix(2,
Polynomial(Integer))) to type Matrix(Matrix(Polynomial(Integer)))
.

Apparently, with generalized matrices interpeter decided that
Matrix(Matrix(Polynomial(Integer))) would be good type. But
'setelt!' can not change type and interpeter has limited
capablitiy for backtracking, so instead of trying someting else
it produces error message.

Note: interpreter bugs are not new, but they are more visible
with new types.

Anyway, currently IMO we can not commit this code, we need solution
to interpreter problems first. We probably could commit change
to MatrixCategory (which is large part of changes), but this is
of limited use.

--
Waldek Hebisch

[Tim Daly]

I can highly recommend creating a "tensor" package (see

https://www.tensorflow.org/)

This has four advantages.

Tensors are much more general than the matrix approach.

The tensorflow package does just about anything and is an extremely flexible

data structure. (see https://www.tensorflow.org/api_docs/python/tf)

If you scroll down to the "Classes" section you can see how this would

easily fit into the category/domain hierarchy.

Tensors also form the basis for most of the GPT work. This might be of great

benefit for future development.

Tensors work well with GPUs, supporting more future development.

The API and Class structures are already designed and widely used.

Tim

[Ralf Hemmecke]

On 07.05.23 04:49, Waldek Hebisch wrote:
> In mail from 4 Nov 2019 Johannes Grabmeier proposed extending
> two dimensional arrays to work as matrices of matrices.

Link:
https://www.mail-archive.com/fricas...@googlegroups.com/msg13038.html

> I have looked at possible ways to obtain desired functionality.

I must admit that I did not completely understand what the expected
result should be. So my guess is that FriCAS should be able to work with
block matrices.

We already have Matrix. That allows any rectangular form of a matrix.
So from a naive point of view matrices can be entries of a "meta" matrix.

Of course, one wants to restrict the shapes of the entries so that the
result can be seen just as one matrix with a substructure.
And certainly, one can have matrices of matrices of matrices.

I actually see no big problem in implementing this, it is just that one
would have to take care of the compatibility of the + and * operations.
The only thing is that this new type constructor shouldn't be (as a
first approach) of type MatrixCategory, since this requires its entry
type to be an abelian monoid. If we do not require this, but only
require an operation +: (R, R) -> R that would be sufficient.
AFAIR Bronstein introduced a so-called ArithmeticType into the
Aldor-Library. that would be a + with no properties. No commutativity or
associativity. Also a bit quesionable, but an idea. For matrices, we
would even run into the problem that + is not defined for all input.
But that is also fine. The operation just fails for block matrices if
the dimensions of the blocks are not appropriate.

> I have now mostly implemented the first possiblity, that is
> generalization of Matrix. AFAICS generalized domain works
> fine. However, it uncovered troubles with interpeter and
> Spad compiler. Namely, Spad compiler has trouble with
> condition like
>
> if UP has zero? : UP -> Boolean then
>
> where UP is domain parameter which satisfies
> UnivariatePolynomialCategory(F) and F is UniqueFactorizationDomain.
> This imples that conditon above is satisfies, but Spad
> compiler can not infer this, so I had to add extra
> conditionals as workarounds.

Oh yes, but that is only a weakness of the compiler and not the main
point why Matrix(Matrix Integer) shouldn´t be possible.

> Interpreter problems are more troubling. First, after generalization
> interpreter can generate Matrix(Symbol), which alone is positive
> fact.

Hmmmm... then you see this different from me. Already in my first
reaction to Johannes' mail I spoke about data domains and arithmetic
domains. When I hear Matrix, I immediately think of being able to add or
multiply these structure. So a two dimensional array of Symbol shold be
ok, but I would require the input type of Matrix to at least provide +
and *.

But OK, if you are in favour of generalizing that requirement, I am not
completely against. Already no we require only AbelianMonoid, but no
multiplicative structure.

> Apparently, with generalized matrices interpeter decided that
> Matrix(Matrix(Polynomial(Integer))) would be good type. But
> 'setelt!' can not change type and interpeter has limited
> capablitiy for backtracking, so instead of trying someting else
> it produces error message.

Admittedly, an error message is not so nice, but as long as a user can
give the right types to the entries and then things work fine, I see no
problem.

In fact, I actually do not understand "setelt! cannot change type". What
exactly is the problem.

> Anyway, currently IMO we can not commit this code, we need solution
> to interpreter problems first.

Can you push the code to a new branch of *your* github fricas repo under
https://github.com/hebisch ?

Thanks
Ralf

[Waldek Hebisch]

On Sat, May 13, 2023 at 04:37:08PM +0200, Ralf Hemmecke wrote:
> On 07.05.23 04:49, Waldek Hebisch wrote:
> > In mail from 4 Nov 2019 Johannes Grabmeier proposed extending
> > two dimensional arrays to work as matrices of matrices.
>
> Link:
> https://www.mail-archive.com/fricas...@googlegroups.com/msg13038.html
>
> > I have looked at possible ways to obtain desired functionality.
>
> I must admit that I did not completely understand what the expected result
> should be. So my guess is that FriCAS should be able to work with block
> matrices.

Well, finding out exact specification for new functionality is
important part of implementation. Clearly, we want to work with
block matrices. But I want various parts to be as general as
possible. And there is question of interaction/cooperation
with existing features.

> We already have Matrix. That allows any rectangular form of a matrix.
> So from a naive point of view matrices can be entries of a "meta" matrix.
>
> Of course, one wants to restrict the shapes of the entries so that the
> result can be seen just as one matrix with a substructure.

That is not so clear: for addition it is enough to have compatible
shapes.

> And certainly, one can have matrices of matrices of matrices.
>
> I actually see no big problem in implementing this, it is just that one
> would have to take care of the compatibility of the + and * operations.
> The only thing is that this new type constructor shouldn't be (as a first
> approach) of type MatrixCategory, since this requires its entry type to be
> an abelian monoid.

Of course when we generalize matrices (as opposed to creating something
different) first step is to generalize MatrixCategory.

> If we do not require this, but only require an operation
> +: (R, R) -> R that would be sufficient.
> AFAIR Bronstein introduced a so-called ArithmeticType into the
> Aldor-Library. that would be a + with no properties. No commutativity or
> associativity. Also a bit quesionable, but an idea.

Well,

with "+" : (R, R) -> R

should be resonably good name for such type, but new word probably
is easier for compiler/interpreter.

> For matrices, we would
> even run into the problem that + is not defined for all input.
> But that is also fine. The operation just fails for block matrices if the
> dimensions of the blocks are not appropriate.
>
> > I have now mostly implemented the first possiblity, that is
> > generalization of Matrix. AFAICS generalized domain works
> > fine. However, it uncovered troubles with interpeter and
> > Spad compiler. Namely, Spad compiler has trouble with
> > condition like
> >
> > if UP has zero? : UP -> Boolean then
> >
> > where UP is domain parameter which satisfies
> > UnivariatePolynomialCategory(F) and F is UniqueFactorizationDomain.
> > This imples that conditon above is satisfies, but Spad
> > compiler can not infer this, so I had to add extra
> > conditionals as workarounds.
>
> Oh yes, but that is only a weakness of the compiler and not the main point
> why Matrix(Matrix Integer) shouldn´t be possible.

Yes.

> > Interpreter problems are more troubling. First, after generalization
> > interpreter can generate Matrix(Symbol), which alone is positive
> > fact.
>
> Hmmmm... then you see this different from me. Already in my first reaction
> to Johannes' mail I spoke about data domains and arithmetic domains. When I
> hear Matrix, I immediately think of being able to add or multiply these
> structure. So a two dimensional array of Symbol shold be ok, but I would
> require the input type of Matrix to at least provide + and *.

Well, in the past you wrote "mathematical domains". I can imagine
mathematical operations on Matrix(Symbol). And it would be natural
for interpreter to convert Matrix(Symbol) to matrix of polynomials.

> But OK, if you are in favour of generalizing that requirement, I am not
> completely against. Already no we require only AbelianMonoid, but no
> multiplicative structure.
>
> > Apparently, with generalized matrices interpeter decided that
> > Matrix(Matrix(Polynomial(Integer))) would be good type. But
> > 'setelt!' can not change type and interpeter has limited
> > capablitiy for backtracking, so instead of trying someting else
> > it produces error message.
>
> Admittedly, an error message is not so nice, but as long as a user can give
> the right types to the entries and then things work fine, I see no problem.
>
> In fact, I actually do not understand "setelt! cannot change type". What
> exactly is the problem.

Interpreter considers 'setelt!' to be two argument operation. Normally
interpreter will coerce arguments to required types and function
selection assumes that if argument types are coercible to needed
type, then they are OK. But 'setelt!' is special, it works by
side effect, so it must work on argument "as given" and can
not coerce it to new type (which would create copy).

> > Anyway, currently IMO we can not commit this code, we need solution
> > to interpreter problems first.
>
> Can you push the code to a new branch of *your* github fricas repo under
> https://github.com/hebisch ?

I have put the code (diff) at

http://www.math.uni.wroc.pl/~hebisch/fricas/matmat.diff

--
Waldek Hebisch

[Ralf Hemmecke]

> Well, finding out exact specification for new functionality is
> important part of implementation. Clearly, we want to work with
> block matrices. But I want various parts to be as general as
> possible. And there is question of interaction/cooperation
> with existing features.

In that case one could even allow block matrices of the following form

a11 :+ matrix [[1]]
a12 := matrix [[2],[3]]
a21 := matrix [[22,33,44]]
ass := matrix [[7]]
mat := matrix [[a11,a22],[a21,a22]]

That doesn´t look like what on usually considers a block matrix, but mat
would be an object that is perfectly fine to do mat + mat or 3*mat.
Multiplication mat*mat should fail due to incompatibilities as would
addiion of incompatible sizes. So we could introduce +: (%, %) -> % as a
partial function that gives an aborts the computation if the arguemnts
are incompatible.

I see this as good and bad.

I usually expect a function to be total, however, already / is a clear
exception and turning the signature of / into
(%, %) -> Union(%,"failed")
is clearly an overkill.

I would in this discussion rather restrict to Matrix and clearly define
what the mathematical structure is. For example, if arguments are
compatible, then addition is allowed and yields a result, if furthermore
the argument type of Matrix(X) is commutative/associative then so is the
addition of matrices. Etc, etc.

I have the impression that you have no problem going this way.

> Well,
>
> with "+" : (R, R) -> R
>
> should be resonably good name for such type, but new word probably
> is easier for compiler/interpreter.

Hmm... you probably started this thread not because the implementation
is problematic, but the interaction with other parts of FriCAS (in
particular the interpreter) is.

> Well, in the past you wrote "mathematical domains". I can imagine
> mathematical operations on Matrix(Symbol). And it would be natural
> for interpreter to convert Matrix(Symbol) to matrix of polynomials.

For me coercion counts as embedding and that should never cause problems
if we even require that to be a homomorphism, but the interpreter also
applies convert and retract (as far as I understand) and that might not
give what the user wanted.

Clearly the type system is good for programming, but I admit that it is
a quite difficult task to guess the types that the user intended. That
the interpreter fails is no surprise. In fact, sometimes I would have
liked more control over the exposed domains.

Although probably unnatural for most users, I would be happy with a
session starting with zero knowledge where I would have to say "import
from Integer" before I can do 1+1.

Side remark: Yes, I would rather like "import from XXX" instead of
)expose XXX.

>> In fact, I actually do not understand "setelt! cannot change type". What
>> exactly is the problem.
>
> Interpreter considers 'setelt!' to be two argument operation. Normally
> interpreter will coerce arguments to required types and function
> selection assumes that if argument types are coercible to needed
> type, then they are OK. But 'setelt!' is special, it works by
> side effect, so it must work on argument "as given" and can
> not coerce it to new type (which would create copy).

Yes. Clear. Do you mean the problem is that setelt! is currently not
treated special? Or what exactly is the problem. Can you give an example?

Ralf

[Waldek Hebisch]

On Sun, May 14, 2023 at 12:26:20AM +0200, Ralf Hemmecke wrote:
> > Well, finding out exact specification for new functionality is
> > important part of implementation. Clearly, we want to work with
> > block matrices. But I want various parts to be as general as
> > possible. And there is question of interaction/cooperation
> > with existing features.
>
> In that case one could even allow block matrices of the following form
>
> a11 :+ matrix [[1]]
> a12 := matrix [[2],[3]]
> a21 := matrix [[22,33,44]]
> ass := matrix [[7]]
> mat := matrix [[a11,a22],[a21,a22]]
>
> That doesn´t look like what on usually considers a block matrix, but mat
> would be an object that is perfectly fine to do mat + mat or 3*mat.

With the patch I gave (after fixing typos) one can create 'mat' as
above. 'mat + mat' works. '3*mat' gives weird error, message is
rather confusing, but reason is simple: Matrix(Matrix(R)) has
multiplication by Matrix(R), but not by R...

> Multiplication mat*mat should fail due to incompatibilities as would addiion
> of incompatible sizes.

Yes.

> So we could introduce +: (%, %) -> % as a partial
> function that gives an aborts the computation if the arguemnts are
> incompatible.

Well, this already is the case, nothing really new here.

> I see this as good and bad.
>
> I usually expect a function to be total, however, already / is a clear
> exception and turning the signature of / into
> (%, %) -> Union(%,"failed")
> is clearly an overkill.
>
> I would in this discussion rather restrict to Matrix and clearly define what
> the mathematical structure is. For example, if arguments are compatible,
> then addition is allowed and yields a result, if furthermore the argument
> type of Matrix(X) is commutative/associative then so is the addition of
> matrices. Etc, etc.

Yes.

> I have the impression that you have no problem going this way.
>
> > Well,
> >
> > with "+" : (R, R) -> R
> >
> > should be resonably good name for such type, but new word probably
> > is easier for compiler/interpreter.
>
> Hmm... you probably started this thread not because the implementation is
> problematic, but the interaction with other parts of FriCAS (in particular
> the interpreter) is.

My comment was about ArithmeticType.

> > Well, in the past you wrote "mathematical domains". I can imagine
> > mathematical operations on Matrix(Symbol). And it would be natural
> > for interpreter to convert Matrix(Symbol) to matrix of polynomials.
>
> For me coercion counts as embedding and that should never cause problems if
> we even require that to be a homomorphism, but the interpreter also applies
> convert and retract (as far as I understand) and that might not give what
> the user wanted.
>
> Clearly the type system is good for programming, but I admit that it is a
> quite difficult task to guess the types that the user intended. That the
> interpreter fails is no surprise. In fact, sometimes I would have liked more
> control over the exposed domains.
>
> Although probably unnatural for most users, I would be happy with a session
> starting with zero knowledge where I would have to say "import from Integer"
> before I can do 1+1.
>
> Side remark: Yes, I would rather like "import from XXX" instead of )expose
> XXX.

The "best" which you can do ATM is

)set expose drop group basic

After that, given

vector([1, 2])

interpreter complains that it can not find exposed 'elt'. But

1 + 1

and

[1, 2]

still work, presumably due to hardcoded info.

> > > In fact, I actually do not understand "setelt! cannot change type". What
> > > exactly is the problem.
> >
> > Interpreter considers 'setelt!' to be two argument operation. Normally
> > interpreter will coerce arguments to required types and function
> > selection assumes that if argument types are coercible to needed
> > type, then they are OK. But 'setelt!' is special, it works by
> > side effect, so it must work on argument "as given" and can
> > not coerce it to new type (which would create copy).
>
> Yes. Clear. Do you mean the problem is that setelt! is currently not treated
> special?

It is not treated special when selecting which operation to use.
Interpeter only realises that it is special when it should coerce
the first argument...

> Or what exactly is the problem. Can you give an example?

Example is from 'fixed.input':

msq := Matrix SquareMatrix(2,POLY INT)
m : msq := zero(2,2)
m(1,1) := matrix([[1,2],[a,b]])

--
Waldek Hebisch