Matrix Symbol

[Matthew Rocklin]

Hi Everyone, 

I'd like to build a symbolic Matrix object. I'm envisioning a non-commutative Symbol with a shape and some optional attributes like is_symmetric, is_positive_definite, etc.... I imagine that Inverse and Transpose would be objects on top of expressions built up of these. 

To be clear this Matrix won't hold any data, it won't be mutable in any way, it'll just be a mathematical symbol. 

It may be that this work conflicts with Sherjil's work on revamping matrices. If that's the case I'll just hold off until some more ideal future when things are quiet. I'm hoping that this is relatively independent from his work.

I was looking through the issues and it appears that similar thoughts have come up in the past. What's the general community wisdom on this topic?

Best, 

-Matt

[Aaron Meurer]

This has already been considered some. See
http://code.google.com/p/sympy/issues/detail?id=1860 for example.

I kind of doubt this would conflict with Sherjil's work, since what
you are suggesting and actual matrices would be quite different, but
Sherjil would have to confirm.

I think if you implement non-commutative symbols with shape (see the
issue I linked to above), that would cover about 90% of use cases. If
you want to symbolically represent more than just matrix
multiplication and addition, you would have to do more than that, but
these would not be difficult to do. For example, A.T would just swap
the shape of A. The inverse would just be handled by Pow, i.e.,
A**-1. I suppose a little more work would need to be done if you want
to support matrices being non-invertible, and that being handled
correctly. But even that would not be hard (just have an attribute .
invertible, which would probably default to True if the shape is
square, and make ._eval_pow, not work for negative powers if it is
False).

The real problem (and I think it's what the main issue is in issue
1860), is how to make Mul handle the shape attribute correctly. The
way it is now, you have to modify Mul.flatten to make it work (also
see issue 1941).

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.
>

[Simon]

There was some discussion of this in the sage group. The code will also work in SymPy:

http://ask.sagemath.org/question/505/symbolic-matrices

(also see http://stackoverflow.com/questions/5708208/symbolic-matrices-in-mathematica-with-unknown-dimensions/5710838)

This type of approach seems to be what you were talking about, and can easily be extended to contain the is_positive_definite type properties.

It needs some work to make inverse work nicely with a square product of rectangular matrices, etc...

[Matthew Rocklin]

@Simon: Thanks for the pointer to those pages. They're good to look through

@Aaron: Hrm, yes getting Muls and Adds and such to recognize matrices
is quite a hurdle. Should I shy away from adding to flatten()?
Alternatively I could subclass all of these into MatrixExpr, MatrixMul
(MatMul?), MatrixAdd, etc... and build in a couple of checks into the
__add__, __mul__ operators before passing up to the superclass method.

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.

> To view this discussion on the web visit
> https://groups.google.com/d/msg/sympy/-/OOTzZmk0RPoJ.

[SherjilOzair]

Hello Matthew,

As Aaron says correctly, work on matrices and work on the matrix symbol is quite separate.

The Matrix class provides fundamental operations like addition, multiplication. We only have to make Add, Mul call the required routines in the Matrix class.

I'm not sure what .flatten does, and my knowledge about the core is only minimal. But I think the only time there is a difference between a regular symbol and a matrix symbol is when expr.subs(x, matrix) is done. Possibly, this could be achieved only by adding some code to the subs function. Instead of focussing on things like is_symmetric, is_positive_definite (which I agree, are useful), I think making expr.subs(x, matrix) work should be the first step. x here could be a regular non-commutative Symbol.

Even if its inelegant, i.e. involves making Add, Mul Matrix-aware, I think it would be good if we could try out stuff making .subs work for matrices by adding code to the present core (.flatten, .subs, etc.). If we succeed, then more code could be added to get more functionality, even if we don't push in stuff to the master. Later if this successful, all the Matrix expressions code could be ported to a separate MatrixMul, MatrixAdd, etc. Making new classes that belong to the core for a small functionality feels like overkill right now.

Have you tried anything out regarding this, Matthew ?

Regards,

Sherjil Ozair

[Matthew Rocklin]

Regarding expr.subs(x,matrix):

This is tricky because Matrices aren't sympifiable. I think you would need to substitute all matrix symbols for matrices at once and then turn the expression into a Matrix object using the Matrix.__add__ etc methods. Or maybe any substituted Matrix would become some sort of ImmutableMatrix object so that it could live inside an Expr.

Regarding other differences:

Being able to substitute computational elements like Matrices/SparseMatrices/LinearOperators into matrix expressions is definitely the main motivating application behind symbolic matrices. This is not my major focus right now though. I'm focusing on turning statistical statements into computational ones. I've been generating integrals for a while now and would like to start generating Matrix Expressions like what's found here. Being able to check the shape of an expression is likely going to be useful for me. I suspect that there will be more features of matrices that I'd like exposed other than just that they're not commutative. Knowing if it is invertible or not is a good example.  

Of course, after I generate and simplify these expressions I'll need to substitute in actual matrices of some form or another, I just haven't gotten there yet :)

Haven't tried anything yet.

-Matt

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To view this discussion on the web visit

> https://groups.google.com/d/msg/sympy/-/uyTjPpzKQhQJ.

[Aaron Meurer]

Unless you can make matrices sympifyable, you probably won't be able
to use .subs. That isn't to say that you can't write your own
specialized function to do the substitution.

There's also some cool stuff you could do here, like implementing a
matrix chain multiplication algorithm
(http://en.wikipedia.org/wiki/Matrix_chain_multiplication) to compute
the most efficient way to multiply a product of matrices ABCD...

Aaron Meurer

[Matthew Rocklin]

Yeah, definitely. I must confess that a hidden passion of mine is optimizing linear algebra (we all have our quirks). I was just looking at Theano a minute ago actually - I think it would be cool to easily dump Matrix expressions onto them.

How should matrix expressions be represented in SymPy though? The way I see it there are three options

1. Leave it as a SymPy Expr, forget shape, transpose, rank, etc... for now. Maybe future SymPy elegance will make clever things possible (such as issue 1941)
2. Enhance SymPy Expr's to play nice with Matrices
3. Subclass a family of MatrixExpr classes that live outside the core

Probably there are things I'm missing but this is how I separate things. Because I'd like this done sooner rather than later I'm obviously in favor of 2 or 3 with a preference for 3. I don't know enough about SymPy however to tell whether this is the "right way" and I'd rather not work on something unless it has a chance at getting in. 

I'll push again for three by saying that there is a lot going on in Matrix Expressions other than just non-commutativity and shape. Inverses, transposes, rank, symmetry, positive_definiteness, conditioning, etc... all play a big role in computational decisions made on matrices. Additionally matrix expressions are ubiquitous and important in the scientific computing world. From my (strongly biased) perspective they deserve special treatment. 

But for my Summer project I just need something. Even 1 might work for me. Does anyone have any thoughts or suggestions?

Best,

-Matt 

[Aaron Meurer]

Option 3 might work. You should talk to Brian Granger about the
advantages and disadvantages of defining your own Expr objects, as
that's what is done in the quantum code.

Note that you should be able to handle the other stuff with the
assumptions (in theory anyway; assumptions are kind of a mess right
now, so I couldn't say how well it would work).

But even if we fix the core and option 1 becomes feasible, the code
can be refactored, so if you want something now, I am +1 to trying
option 3. I would rather prefer to avoid option 2, unless you are
going to fix it in a more general way (like issue 1941).

Aaron Meurer

[Brian Granger]

Matthew,

I *strongly* recommend against subclassing Mul/Pow/Add as you will
find that it is *impossible* even with op_priority to have your
Mul/Pow/Add used all of the time. In a variety of contexts, you will
*want* your subclasses to be used, but sympy's base Mul/Add/Pow will
be used. There is no way to prevent this with the current design.
Just use the sympy's existing Mul/Add/Pow and live with its
limitations for now.

Cheers,

Brian

On Tue, Jun 28, 2011 at 9:00 AM, Matthew Rocklin <mroc...@gmail.com> wrote:
> There is a barebones implementation of a Matrix Expression class here:
> https://github.com/mrocklin/sympy/tree/matrix_expr/sympy/matrices
> Comments/criticism welcome.

--
Brian E. Granger
Cal Poly State University, San Luis Obispo
bgra...@calpoly.edu and elli...@gmail.com

[Vinzent Steinberg]

On Jun 27, 5:54 pm, Matthew Rocklin <mrock...@gmail.com> wrote:
> Regarding expr.subs(x,matrix):
>
> This is tricky because Matrices aren't sympifiable. I think you would need
> to substitute all matrix symbols for matrices at once and then turn the
> expression into a Matrix object using the Matrix.__add__ etc methods. Or
> maybe any substituted Matrix would become some sort of ImmutableMatrix
> object so that it could live inside an Expr.

I would rather use lambdify than subs. Actually, I think this might
work already.

Vinzent

[Vinzent Steinberg]

On Jun 28, 4:32 am, Matthew Rocklin <mrock...@gmail.com> wrote:
> Yeah, definitely. I must confess that a hidden passion of mine is optimizing
> linear algebra (we all have our quirks). I was just looking at Theano a
> minute ago actually - I think it would be cool to easily dump Matrix
> expressions onto them.
>
> How should matrix expressions be represented in SymPy though? The way I see
> it there are three options
>

>    1. Leave it as a SymPy Expr, forget shape, transpose, rank, etc... for

>    now. Maybe future SymPy elegance will make clever things possible (such as
>    issue 1941)

>    2. Enhance SymPy Expr's to play nice with Matrices
>    3. Subclass a family of MatrixExpr classes that live outside the core

>
> Probably there are things I'm missing but this is how I separate things.
> Because I'd like this done sooner rather than later I'm obviously in favor
> of 2 or 3 with a preference for 3. I don't know enough about SymPy however
> to tell whether this is the "right way" and I'd rather not work on something
> unless it has a chance at getting in.
>
> I'll push again for three by saying that there is a lot going on in Matrix
> Expressions other than just non-commutativity and shape. Inverses,
> transposes, rank, symmetry, positive_definiteness, conditioning, etc... all
> play a big role in computational decisions made on matrices. Additionally
> matrix expressions are ubiquitous and important in the scientific computing
> world. From my (strongly biased) perspective they deserve special
> treatment.

I think all this '.is_*' stuff should rather be implemented using the
new assumption system (not sure about shape, maybe this should really
go into the core). If we use noncommutative symbols, I think we can
avoid messing around with Mul and friends.

A.T could be implemented as a unary function.

Vinzent

[Matthew Rocklin]

There is a barebones implementation of a Matrix Expression class here:

https://github.com/mrocklin/sympy/tree/matrix_expr/sympy/matrices

Comments/criticism welcome. 

[Matthew Rocklin]

@Brian - Thanks for the heads up Brian. I'll see what I can do with option (1). My short term solution was to start a "matrixify" function a la sympify. It would probably be too annoying to use everywhere though. 

@Vinzent - Where is a good place to start learning about the new assumption system (or the old one... I'm not up to speed on these)

--

[Amit]

Hi,

I am not familiar with the internals of sympy. But I suggest that if
you start working on the implementation of symbolic matrices, you
should take into consideration more complicated operators like
differentiation.
'The Matrix Cookbook' has many matrix equalities that maybe can be
implemented using some kind of pattern recognition.

Amit

[Aaron Meurer]

As I pointed out in the other thread, non-commutative differentiation
already works in SymPy, so doing this should not be difficult.

Aaron Meurer

[Alan Bromborsky]

Differentiation would only work with a scalar (communicative)
differentiation operator. If the matrix function is a function of a
vector or matrix one would have to define the directional derivative for
each case (which would be a scalar differential operator) and use the
results of that operation to determine the properties of a vector or
matrix derivative. Note that determining the operator properties would
also require a definition for the scalar product of vectors and matrices.

Consider the vector directional derivative of a matrix M that is the
function of a vector v, M(v), then if a is an arbitrary vector (LaTeX
expression) the definition of the directional derivative would be

a\cdot\nabla M(v) \equiv \lim_{h \rightarrow 0}\frac{M(v+ha)-M(v)}{h}

from this the properties of \nabla M(v) could be determined and if v is
expanded in a arbitrary basis the \nabla operator could also be
expanded. A similar treatment is possible for a matrix that is a
function of a matrix if the scalar product of two matrices is defined.

[Matthew Rocklin]

I agree that support for derivatives on matrices is important; I would like this myself. I haven't thought much about it in the context of SymPy before though so thank you for bringing it up.

I haven't solidified my understanding of this problem but it seems like there are a few concepts of a derivative here. 

• Given a matrix expression we can differentiate with respect to the various matrices themselves. This is likely very similar to Aaron's example using stock SymPy with non-commutative symbols. This should (I think) work out of the box
• Given a function on vector valued inputs with possible vector valued outputs we can define directional derivatives. I think this is what Alan is talking about. In this situation the objects can easily become high rank (Hessians of vector-valued functions are not matrices). This leads me to the idea that we should consider mathematical Tensors rather than matrices. 

The tensor vs matrix issue makes me nervous. There is a lot of special stuff happening on rank two tensors (matrices) that we rarely care about at higher ranks. Maybe this work should be done on a Tensor class and Matrices should subclass Tensors? This is getting beyond the scope of what I was looking for with a Matrix Symbol. I probably won't pursue this enhancement in the near future but would be happy to support someone else's effort. 

For the moment I'm not working on Matrix Expressions actually. I'm a bit stuck on how to proceed and would welcome any suggestions. The best idea I have now is to insert symbols into standard sympy Exprs and have aspects like shape and is_invertible be functions which are called on the Expr tree rather than fields or methods of the object. This will fail to raise exceptions when illegal operations are performed but should get the job done. The Indexed class is somewhat similar in flavor. 

On Thu, Jun 30, 2011 at 7:05 AM, Alan Bromborsky <abr...@verizon.net> wrote:

Differentiation would only work with a scalar (communicative) differentiation operator.  If the matrix function is a function of a vector or matrix one would have to define the directional derivative for each case (which would be a scalar differential operator) and use the results of that operation to determine the properties of a vector or matrix derivative.  Note that determining the operator properties would also require a definition for the scalar product of vectors and matrices.

Consider the vector directional derivative of a matrix M that is the function of a vector v, M(v), then if a is an arbitrary vector (LaTeX expression) the definition of the directional derivative would be

a\cdot\nabla M(v) \equiv  \lim_{h \rightarrow 0}\frac{M(v+ha)-M(v)}{h}


from this the properties of \nabla M(v) could be determined and if v is expanded in a arbitrary basis the \nabla operator could also be expanded.  A similar treatment is possible for a matrix that is a function of a matrix if the scalar product of two matrices is defined.

On 06/30/2011 04:20 AM, Aaron Meurer wrote:

As I pointed out in the other thread, non-commutative differentiation
already works in SymPy, so doing this should not be difficult.

Aaron Meurer

sympy+unsubscribe@googlegroups.com.

For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.

To unsubscribe from this group, send email to sympy+unsubscribe@googlegroups.com.

For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.

To unsubscribe from this group, send email to sympy+unsubscribe@googlegroups.com.

[Aaron Meurer]

This is something that you would (hopefully) be able to do with the
new assumptions system. In other words, without having to modify the
core, you should be able to say ask(Q.invertable(A*B)) and it would
determine it based on whether you set A and B to be invertible. Shape
should work too, though I'm not sure if the system can currently
handle non-boolean assumptions (someone else will have to fill in
here).

By the way, is_invertible is perhaps something that could be
implemented in the core directly, so you could have support for
non-invertible symbols in any case. It should just be a matter of
making ._eval_power do the right thing in any case.

Aaron Meurer

>
>
> On Thu, Jun 30, 2011 at 7:05 AM, Alan Bromborsky <abr...@verizon.net>
> wrote:
>>
>> Differentiation would only work with a scalar (communicative)
>> differentiation operator.  If the matrix function is a function of a vector
>> or matrix one would have to define the directional derivative for each case
>> (which would be a scalar differential operator) and use the results of that
>> operation to determine the properties of a vector or matrix derivative.
>>  Note that determining the operator properties would also require a
>> definition for the scalar product of vectors and matrices.
>>
>> Consider the vector directional derivative of a matrix M that is the
>> function of a vector v, M(v), then if a is an arbitrary vector (LaTeX
>> expression) the definition of the directional derivative would be
>>
>> a\cdot\nabla M(v) \equiv  \lim_{h \rightarrow 0}\frac{M(v+ha)-M(v)}{h}
>>
>>
>> from this the properties of \nabla M(v) could be determined and if v is
>> expanded in a arbitrary basis the \nabla operator could also be expanded.  A
>> similar treatment is possible for a matrix that is a function of a matrix if
>> the scalar product of two matrices is defined.
>>
>>
>>
>> On 06/30/2011 04:20 AM, Aaron Meurer wrote:
>>>
>>> As I pointed out in the other thread, non-commutative differentiation
>>> already works in SymPy, so doing this should not be difficult.
>>>
>>> Aaron Meurer
>>>

>>>>>> sympy+un...@googlegroups.com.

>>>>>> For more options, visit this group at
>>>>>> http://groups.google.com/group/sympy?hl=en.
>>>>
>>>> --
>>>> You received this message because you are subscribed to the Google
>>>> Groups "sympy" group.
>>>> To post to this group, send email to sy...@googlegroups.com.
>>>> To unsubscribe from this group, send email to

>>>> sympy+un...@googlegroups.com.

>>>> For more options, visit this group at
>>>> http://groups.google.com/group/sympy?hl=en.
>>>>
>>>>
>>
>> --
>> You received this message because you are subscribed to the Google Groups
>> "sympy" group.
>> To post to this group, send email to sy...@googlegroups.com.
>> To unsubscribe from this group, send email to

>> sympy+un...@googlegroups.com.

>> For more options, visit this group at
>> http://groups.google.com/group/sympy?hl=en.
>>
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to

> sympy+un...@googlegroups.com.

[Matthew Rocklin]

Where is the best place to read about the new assumptions system?

[Ronan Lamy]

Le jeudi 30 juin 2011 à 13:25 -0500, Matthew Rocklin a écrit :
> Where is the best place to read about the new assumptions system?

I'm afraid that the best place is the source code in sympy/assumptions/.
I'm not aware of any comprehensive and current write-up on the new
assumptions.

> >>> On Thu, Jun 30, 2011 at 1:58 AM, Amit<amit...@gmail.com>

> +unsub...@googlegroups.com.

> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.
>
>
>
>
>
> --
> You received this message because you are subscribed to the Google
> Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to sympy

> +unsub...@googlegroups.com.

[Aaron Meurer]

There's some info at http://docs.sympy.org/0.7.0/modules/assumptions.html.

Aaron Meurer

[Matthew Rocklin]

It looks like the assumptions system would be a good place for the fancy matrix questions (positive definite, invertible, etc...). I'm not sure that it's appropriate for shape though. 

I'm going to boldly propose that we add shape as a property to the core Expr class. It won't be checked as operations happen (to optimize speed for the common non-matrix case) but will be computed on demand and raise errors at that time if things are misaligned. I might also push for a is_Matrix flag. 

i.e.

>>> A = MatrixSymbol('A', 3, 4)

>>> B = MatrixSymbol('B', 4, 7)

>>> C = 2*A*B # Nothing matrix-y happens here

>>> C.shape # Shape is checked here

(3, 7)

>>> D = 1 + B*A # This passes even though it's invalid

>>> D.shape # Error is found when asking for the shape

ShapeError(....

>>> print Symbol('x').shape # Normal expressions have no shape

None

I don't think this will affect computation at all for non-matrix situations. It adds a minimal amount of code complexity to the Add, Mul, Pow, etc... classes. 

I'm going to write something up for review later tonight. Please someone stop me if this is a poor idea. 

>>         >>> On Thu, Jun 30, 2011 at 1:58 AM, Amit<amiti.bo@gmail.com>

[Brian Granger]

On Sun, Jul 3, 2011 at 1:41 PM, Matthew Rocklin <mroc...@gmail.com> wrote:
> It looks like the assumptions system would be a good place for the fancy
> matrix questions (positive definite, invertible, etc...). I'm not sure that
> it's appropriate for shape though.
> I'm going to boldly propose that we add shape as a property to the core Expr
> class. It won't be checked as operations happen (to optimize speed for the
> common non-matrix case) but will be computed on demand and raise errors at
> that time if things are misaligned. I might also push for a is_Matrix flag.
> i.e.

I am -1 on putting shape into Expr itself. It is too specialized to
be there. It should only be in the classes that actually have shape
information. I agree though that we need to enhance sympy to know
about such things, but that logic type of logic should be in Mul and
should not be narrowly focused on shape but should include the
handling of all types that need custom logic upon multiplication.

As I have said before, I would forget about shape checking *entirely*
at this point in the game though so you can focus on the more
interesting and relevant aspects of this work. Once sympy is extended
in the right way, you will get shape checking almost for free...

Cheers,

Brian

>> >>         >>> On Thu, Jun 30, 2011 at 1:58 AM, Amit<amit...@gmail.com>

--

[Matthew Rocklin]

Hrm, I'm feeling blocked on this issue. I'd really like to be able to generate matrix expressions for the next section of my project and I'd like them to be able to do some things beyond what should go into the standard Expr class. 

Do I have any options here other than to wait until the broader SymPy community decides what it's going to do about smarter handling of subclassed Expr's? 

Brian, you strongly recommend against subclassing Expr. Have you come up with a decent temporary alternative for the physics branch?

I could look into fixing the larger SymPy issue. Is this likely beyond my scope?

If anyone has suggestions I'm happy to hear them. 

Right now I think my best course of action is to continue subclassing Expr very very carefully. 

-Matt

>> >>         >>> On Thu, Jun 30, 2011 at 1:58 AM, Amit<amiti.bo@gmail.com>

[Brian Granger]

Matthew,

On Thu, Jul 7, 2011 at 10:18 AM, Matthew Rocklin <mroc...@gmail.com> wrote:
> Hrm, I'm feeling blocked on this issue. I'd really like to be able to
> generate matrix expressions for the next section of my project and I'd like
> them to be able to do some things beyond what should go into the standard
> Expr class.

Other than shape checking, what types of things do you need to do?

> Do I have any options here other than to wait until the broader SymPy
> community decides what it's going to do about smarter handling of subclassed
> Expr's?
> Brian, you strongly recommend against subclassing Expr. Have you come up
> with a decent temporary alternative for the physics branch?

I may have said confusing things previously. I do think that
subclassing Expr is OK, just to override methods like
__mul__/__add__/__pow__ to provide custom Mul/Add/Pow logic. Here is
how we handle this...

We simply put the additional logic in standalone functions that walk
an expression tree and perform the needed logic. Something like this:

A = Matrix('A', (4,4))
B = Matrix('B', (4,4))

c = A*B
check_shape(c)

Here Matrix would be a custom subclass of Expr, that has the shape
information in args, but it wouldn't have any custom Mul/Pow/Add logic
implemented.

We have found this approach to be extremely simple to implement and
fully functional as well. For details of how we do this in the
quantum stuff see:

For our subclasses of Expr:

sympy.physics.quantum.qexpr
sympy.physics.quantum.state
sympy.physics.quantum.operator

For our functions that walk expressions:

sympy.physics.quantum.qapply
sympy.physics.quantum.represent

> I could look into fixing the larger SymPy issue. Is this likely beyond my
> scope?

Yes, I think it would be a much larger effort than you project would allow for.

Cheers,

Brian

>> >> >> Amit<amit...@gmail.com>

[Matthew Rocklin]

"Other than shape checking, what types of things do you need to do?"

Good question, one to which I do not have an exhaustive answer. So far Inverses, Transpose, and symbolic block matrices are on the list. I expect this to grow as my use of it continues. This sort of thing is likely to be useful for me outside of my GSoC project and I'd like to think others might find it useful and contribute as well. 

I think it would be good for the growth of symbolic matrix expressions to not be constrained by the purity of the core. 

I do like the idea of standalone functions to check properties inside SymPy Expr's. This is what I do with Random Expressions (my GSoC project). My personal intuition is that Matrix operations are sufficiently distinct to warrant a further deviation from SymPy Expr's though. This is just intuition though.

I suspect that this is a topic for which there exist many diverse opinions. What is the best way to balance making progress on my project and building something that the community will agree with? There are a number of old e-mails and issues about immutable matrices and such. Should I start an issue or new e-mail thread to get input?

-Matt

>> >> >> Amit<amiti.bo@gmail.com>

[Aaron Meurer]

How do you envision symbolic block matrices working? Inversers and
transposes are trivial to implement (see my previous replies to this
thread).

Aaron Meurer

>> >> >> >> Amit<amit...@gmail.com>

[Matthew Rocklin]

In an ideal world they would probably subclass an immutable matrix and use an immutable matrix for storage. This brings up deeper questions. 

My main thesis here though is that we should start a matrix expression section of SymPy because there's lots of good math there and it's tricky to integrate it cleanly into pure SymPy expressions. An experimental MatrixExpr codebase would allow us to develop a bit more freely until SymPy Expr evolves to allow easy integration. I'm just dreaming here though. I don't have a firm picture of how all of this should fit into SymPy and rely on you guys to set me straight.  

Block matrices is just the example that came up yesterday - I'd like to write down a more general version of this. I suspect that there are many more things that I (and others) might want to do with symbolic matrix expressions which are tricky to integrate into pure SymPy expressions. 

>> >> >> >> Amit<amiti.bo@gmail.com>

[Aaron Meurer]

On Thu, Jul 7, 2011 at 8:03 PM, Matthew Rocklin <mroc...@gmail.com> wrote:
> In an ideal world they would probably subclass an immutable matrix and use
> an immutable matrix for storage. This brings up deeper questions.
> My main thesis here though is that we should start a matrix expression
> section of SymPy because there's lots of good math there and it's tricky to
> integrate it cleanly into pure SymPy expressions. An experimental MatrixExpr
> codebase would allow us to develop a bit more freely until SymPy Expr
> evolves to allow easy integration. I'm just dreaming here though. I don't
> have a firm picture of how all of this should fit into SymPy and rely on you
> guys to set me straight.

Well, you're right that it would be good to have an experimental code
base. You will learn much better what should and should not be done
by trying it and seeing what does and does not work than you will
sitting here discussing it. Those of us with more experience can
definitely direct you away from things that we can already see won't
work, but the only way to know the way you *should* do it is to try
implementing it.

Block matrices would be much more work than everything else, because
now you're combining your symbolic matrices with actual matrices,
which themselves would have to be extended to support such things.
Maybe MatrixExpr should just be added as a ground type to Matrix.
That might not be so easy to do, though, since most stuff isn't
written to take anything but shape 1x1 into account. I'd see what
Sherjil thinks about this.

Aaron Meurer

>> >> >> >> >> Amit<amit...@gmail.com>

[Sherjil Ozair]

Implementing symbolic block matrices by using matrices as groundtypes of Matrix is a reasonable idea. But, with the way I have planned groundtypes, it could only support homogenous blocks, i.e. the shape of each block should be the same. I can't think of an elegant way to implement non-homogenous block matrices using the groundtypes method, though.

- Sherjil Ozair

>> >> >> >> >> Amit<amiti.bo@gmail.com>

[SherjilOzair]

There is something I'm doing as part of my project which maybe be
useful. I'm implementing a Matrix_ wrapper class which will wrap over
low-level matrices. Its being written to replace the current Matrix
class. Algorithmic code and user-level code is being separated into
different classes. Currently, we have three internal matrices, the
DenseMatrix( a modified form of the current Matrix class), DOKMatrix
and LILMatrix. These three are essentially Data internal matrices.
The MatrixSymbol can be another object that Matrix_ (later Matrix)
would use as internal.

Of what I understand from this discussion about the Matrix Symbol, we
need to have an object that will be treated like a Matrix everywhere,
but without the internal data. Space is being made in the Matrix
module for such an object, but for it work nicely it needs to interact
nicely with Expr objects. If its a problem subclassing Expr and
friends, I'm +1 about making a separate codebase for matrix
expressions. It is, after all, a different algebra.

I'm not very familiar with Expr, Add, etc. but I think a separate
algebra MatrixExpr, MatrixAdd, etc. could be developed by copy-pasting
some code from Expr, Add, and modifying according to need. I say 'copy-
paste' and not subclassing as the code would be need to be modified
much before it can work on matrices.

Presently, I'm quite busy with this new encapsulation class Matrix_
and can only think of working fulltime on MatrixSymbol after I've
merged in my work on DOKMatrix, LILMatrix and the ongoing work on
Matrix_.

Matthew, you could set up a wiki page on 'matrix expressions' with
some example code in it. 'symbolic matrices' can mean a few number of
things, and is somewhat vague. It would be good if everyone could see
what sort of concrete things are expected from matrix expressions.
Approaches to implement Matrix Algebra, with pros and cons listed,
could be added there for discussion, suggestion and improvement.

Regards,
Sherjil Ozair

[Aaron Meurer]

Yes, that's what I figured. Perhaps rather than being a ground type,
it should be implemented at the level of LIL and DOK (what do you call
this?). In other words, have a special representation of matrices that
works with block matrices. It could probably be dense in nature,
since block matrices aren't generally very big. This could define the
methods that make sense on block matrices (like multiplication if the
block shapes work out), and raise NotImplementedError for the rest (or
however the other representations do this).

Aaron Meurer

>> >> >> >> >> >> Amit<amit...@gmail.com>

[Aaron Meurer]

On Fri, Jul 8, 2011 at 4:00 AM, SherjilOzair <sherji...@gmail.com> wrote:
> There is something I'm doing as part of my project which maybe be
> useful. I'm implementing a Matrix_ wrapper class which will wrap over
> low-level matrices. Its being written to replace the current Matrix
> class. Algorithmic code and user-level code is being separated into
> different classes. Currently, we have three internal matrices, the
> DenseMatrix( a modified form of the current Matrix class), DOKMatrix
> and LILMatrix. These three are essentially Data internal matrices.
> The MatrixSymbol can be another object that Matrix_ (later Matrix)
> would use as internal.
>
> Of what I understand from this discussion about the Matrix Symbol, we
> need to have an object that will be treated like a Matrix everywhere,
> but without the internal data. Space is being made in the Matrix
> module for such an object, but for it work nicely it needs to interact
> nicely with Expr objects. If its a problem subclassing Expr and
> friends, I'm +1 about making a separate codebase for matrix
> expressions. It is, after all, a different algebra.
>
> I'm not very familiar with Expr, Add, etc. but I think a separate
> algebra MatrixExpr, MatrixAdd, etc. could be developed by copy-pasting
> some code from Expr, Add, and modifying according to need. I say 'copy-
> paste' and not subclassing as the code would be need to be modified
> much before it can work on matrices.

As has been discussed before, this is a bad idea. Just subclass Expr.

Aaron Meurer

>
> Presently, I'm quite busy with this new encapsulation class Matrix_
> and can only think of working fulltime on MatrixSymbol after I've
> merged in my work on DOKMatrix, LILMatrix and the ongoing work on
> Matrix_.
>
> Matthew, you could set up a wiki page on 'matrix expressions' with
> some example code in it. 'symbolic matrices' can mean a few number of
> things, and is somewhat vague. It would be good if everyone could see
> what sort of concrete things are expected from matrix expressions.
> Approaches to implement Matrix Algebra, with pros and cons listed,
> could be added there for discussion, suggestion and improvement.
>
> Regards,
> Sherjil Ozair
>

[Matthew Rocklin]

Subclassing Expr has some issues as well. This is what Brian was referring to. Within all of our code we use Add and Mul and don't check if instead we should use some subclass of Add or subclass of Mul. If I feed a matrix expression into these objects then the special matrix structure is lost. This happens if, for example, you call simplify on a matrix expression. 

I think I can get around this though with a few well placed "matrixify" functions. Matrixify is a function which goes through the expression tree and makes appropriate fixes. I've had good success so far with a very very basic Matrixify function. 

Brian, did you have particular horror stories trying to subclass Expr? I'm enthusiastic about my approach but you seemed to have a bad experience. Can you suggest difficult test cases that I should think about?

[Aaron Meurer]

One thing that I have noticed with regard to overriding __mul__ and
__rmul__ (for example) is that you can never completely control what
happens to your class because of association. For example, suppose
that A and B are MatrixExprs and x is some Expr object (say, a
Symbol). Suppose that A*B should give a ShapeError, and that x is a
scalar, which does not matter where it is in the expression. Then you
have to program it so that all of the following give ShapeError:

1. x*(A*B)
2. (x*A)*B
3. (A*x)*B
4. A*(x*B)
5. (A*B)*x
6. A*(B*x)

Let us look at these. Obviously, 1 and 5 will work, since you have
complete control over A.__mul__(B). Similarly, in 4 and 6, you will
end up with A.__mul__(Mul(x, B)), but this is not a problem, since you
can make A.__mul__ check for the left-most noncommutative in a Mul.

But what about 2 and 3? Here, we have __mul__ being called by a Expr
type, namely, Mul(x, A). You want to get B.__rmul__(Mul(x, A)), since
you have implemented shape checking logic in MatrixExpr.__(r)mul__.
But b.__rmul__(a) is called with a*b only if a.__mul__(b) returns
NotImplemented. So, basically, we need Mul.__mul__(MatrixExpr) to
return NotImplemented. Well, let us look at the code for Mul.__mul__
(actually Expr.__mul__):

@_sympifyit('other', NotImplemented)
@call_highest_priority('__rmul__')
def __mul__(self, other):
return Mul(self, other)

Well, this is not very helpful, because the logic is buried in some
decorators. So let us look at the code for those decorators. I will
not bore you by pasting the entire code here (you can see it in
sympy/core/decorators.py), but basically, _sympifyit makes tries to
sympify other, and makes __mul__ return NotImplemented if it fails. I
am not delving into this, because sympify would not fail for
MatrixExpr (otherwise, we would actually have a problem doing any form
of x*A*B).

So call_highest_priority seems to be a better bet. This decorator I
will also not paste here, but it basically lets you define
._op_priority on your object, and if it is greater than Expr's (which
is quite arbitrarily set to 10.0), then it will call
other.__rmul__(self).

But wait, there's more. _op_priority can make the above cases work,
but * is not the only way that SymPy multiplies things. You will also
have to deal with variations on Mul(x, A, B). Mul completely ignores
_op_priority. Unfortunately, even if this may seem like a more
esoteric way to multiply things that you can just recommend users
avoid, it is used internally a lot, because Mul(*args) is more
efficient than reduce(operator.mul, args).

Thus, you see that it is quite impossible to make things work 100% of
the time without modifying the core. And actually, because of the Mul
thing that would not work at all and that is called by so many core
functions and methods, you will not even get something like things
working 90% of the time, but instead things will break when used a
certain way, and you will have hard to track bugs.

Aaron Meurer

[Matthew Rocklin]

I have op_priority set to 11 so that basic arithmetic operations do use MatMuls and MatAdds. MatMul and MatAdd check shape, then call regular Mul/Add __new__ methods to take advantage of the already built flattening functions. After that's done they go through and change classes back to MatMul/MatAdd in case a normal Mul/Add was called within the __new__. So far this has worked very well. All of the associativity cases 1-6 throw the correct error. The more complex expressions that I've tried have also worked fine. 

Calling Mul or Add explicitly does break things. There's no entry point for me to insert clever logic as I don't want to edit core functions. My solution here is simply to say "Use MatMul or MatAdd instead, or, if you must use Mul/Add, call matrixify on the expression afterwards". Matrixify will be guaranteed to return a Matrix expression if MatrixSymbols exist within a normal Expression. I plan to put it in a few key functions like simplify and expand. 

Is this matrixify solution too much of a hack? Hopefully I can hide it away within a few key functions so that the user never needs to know about it unless they explicitly use Add or Mul. If they're using Add or Mul then I expect them to be clever enough to handle matrixify. 

I'm working on this branch . 

[Vinzent Steinberg]

So, if I understood correctly, we would just have to make Mul() raise
a NotImplementedError if it encounters user-defined types?

Vinzent

[Brian Granger]

Matthew,

Arron has provided a fantastic summary of the issues involve. The
summary is that even if you subclass Add/Mul/Pow, sympy will end up
creating expressions that don't use your subclasses everywhere. This
can't be fixed without changing the core, which is outside the scope
of your project. I know I sound like a broken record, but you should
resist the temptation to define your own Add/Mul/Pow classes.

Cheers,

Brian

--

[Brian Granger]

Matthew,

On Wed, Jul 13, 2011 at 6:39 AM, Matthew Rocklin <mroc...@gmail.com> wrote:
> I have op_priority set to 11 so that basic arithmetic operations do use
> MatMuls and MatAdds. MatMul and MatAdd check shape, then call regular
> Mul/Add __new__ methods to take advantage of the already built flattening
> functions. After that's done they go through and change classes back to
> MatMul/MatAdd in case a normal Mul/Add was called within the __new__. So far
> this has worked very well. All of the associativity cases 1-6 throw the
> correct error. The more complex expressions that I've tried have also worked
> fine.
> Calling Mul or Add explicitly does break things. There's no entry point for
> me to insert clever logic as I don't want to edit core functions. My
> solution here is simply to say "Use MatMul or MatAdd instead, or, if you
> must use Mul/Add, call matrixify on the expression afterwards". Matrixify
> will be guaranteed to return a Matrix expression if MatrixSymbols exist
> within a normal Expression. I plan to put it in a few key functions like
> simplify and expand.
>
> Is this matrixify solution too much of a hack? Hopefully I can hide it away
> within a few key functions so that the user never needs to know about it
> unless they explicitly use Add or Mul. If they're using Add or Mul then I
> expect them to be clever enough to handle matrixify.
> I'm working on this branch .

As I understand it Matrixify goes through an expression and makes sure
that your Mul/Add/Pow subclasses are used. If that is the case, I
advise against it because the second a user starts to do things with
the expression they will have to call Matrixify again. Again, I would
just use the default Mul/Add/Pow and provide standalone functions like
check_shape that perform a specific type of action on an expression.

Cheers,

Brian

--

[Ronan Lamy]

Le mercredi 13 juillet 2011 à 10:51 -0700, Brian Granger a écrit :
> Matthew,
>
> Arron has provided a fantastic summary of the issues involve. The
> summary is that even if you subclass Add/Mul/Pow, sympy will end up
> creating expressions that don't use your subclasses everywhere. This
> can't be fixed without changing the core, which is outside the scope
> of your project. I know I sound like a broken record, but you should
> resist the temptation to define your own Add/Mul/Pow classes.

I'm not convinced. What other solution is there?

Also, creating Add/Mul/Pow subclasses looks like the obvious way of
extending sympy, so we want to make it work at some point. Creating
MatrixAdd, etc. will help a lot with that by highlighting all the places
where we assume things like a + b == Add(a, b). So if Matthew thinks he
can make it work without too much pain, compared to the alternatives,
then I think he should go for it.

[Brian Granger]

On Wed, Jul 13, 2011 at 11:19 AM, Ronan Lamy <ronan...@gmail.com> wrote:
> Le mercredi 13 juillet 2011 à 10:51 -0700, Brian Granger a écrit :
>> Matthew,
>>
>> Arron has provided a fantastic summary of the issues involve.  The
>> summary is that even if you subclass Add/Mul/Pow, sympy will end up
>> creating expressions that don't use your subclasses everywhere.  This
>> can't be fixed without changing the core, which is outside the scope
>> of your project.  I know I sound like a broken record, but you should
>> resist the temptation to define your own Add/Mul/Pow classes.
>
> I'm not convinced. What other solution is there?

For now just don't do it. The quantum stuff is equally complicated as
the matrix symbols stuff (in fact from a mathematical perspective,
they are fomally equivalent) and we are surviving just fine without
it. Last summer we started out by subclassing Add/Mul/Pow for the
quantum stuff. For weeks we kept finding things that were broken. It
was horrible and we realized that i) we would end up rewriting much of
sympy and 2) it would still be fundamentally broken. Once we gave up
on that are just used sympy Mul/Add/Pow we started to move 10 times
faster and had very few problems. We did have to give up on our grand
visions of doing custom logic in Mul/Add/Pow, but we have moved that
logic into standalone functions. Those functions are quite simple
though and work just great.

> Also, creating Add/Mul/Pow subclasses looks like the obvious way of
> extending sympy, so we want to make it work at some point. Creating
> MatrixAdd, etc. will help a lot with that by highlighting all the places
> where we assume things like a + b == Add(a, b). So if Matthew thinks he
> can make it work without too much pain, compared to the alternatives,
> then I think he should go for it.

It is not at all obvious to me that the model should be that
Add/Mul/Pow subclasses are encouraged. I think a much better approach
would be closer to what Mathematica has, where objects can define
rules that are applied when the objects are multiplied. IOW,
Mul/Add/Pow themselves should have extension APIs, but those extension
APIs should not be handled through subclassing. But again, this is a
very difficult problem and it is not clear to me what the solution is.
In either case, I stlll feel it is beyond the scope of the current
work.

Cheers,

Brian

[Matthew Rocklin]

What would be really useful here are de-motivational examples. What are some common operations that can't be done with this model? These would make fantastic test cases to see if sub-classing in SymPy is feasible or not. 

Aaron brought up the associativity tests. These work just fine. What are some more complex tasks?

An example that comes to mind is calling simplify on a matrix expression. We have lots of other *simp functions. Adding a matsimp function seems normal. Is simplify the kind of function for which it would be acceptable to add a "if expr.is_Matrix: return matsimp(expr)" statement? If this sort of thing isn't acceptable and if we want to allow the user to call the standard simplify function on matrix expressions then we have a problem. 

[Aaron Meurer]

Well, the dispatch solution only works for op.whatever, but as we
know, reduce(op.whatever, list_of_args) is asymptotically less
efficient than Whatever(*list_of_args).

Actually, half of the problem isn't so much making Add/Mul/etc. work
with dispatch rules, but making it *efficiently* work with dispatch
rules.

Aaron Meurer

[Andy Ray Terrel]

My suggestion would be to continue with the matrixify solution and
write functions to fix up the expression tree as need be. The whole
issue of making Add/Mul/Pow extensible is separate but the logic can
be transferred pretty easily.

-- Andy

[Ondrej Certik]

On Fri, Jul 15, 2011 at 9:17 AM, Andy Ray Terrel <andy....@gmail.com> wrote:
> My suggestion would be to continue with the matrixify solution and
> write functions to fix up the expression tree as need be.  The whole
> issue of making Add/Mul/Pow extensible is separate but the logic can
> be transferred pretty easily.

Exactly. Just get the job done, implement the necessary code, and
whether or not Add/Mul/Pow classes can be extended is a separate
issue, and I would suggest not to, to keep things simple in the core,
so that we can refactor once we finish the assumptions refactoring.

Ondrej

[Mateusz Paprocki]

Hi,

btw. Here is my interpretation of matrix expressions: https://github.com/mattpap/sympy/tree/matrix_expr. Currently it doesn't implement matrix expression arithmetic, but most certainly I'm not going to allow a*M or A + B to "doit" by default, rather I will implement this within expand() framework. What I need are semi evaluated matrix expressions, so det(A^T) can be rewritten to det(A), but a*M must remain unevaluated. This branch is very preliminary, there are no tests for the new stuff and I didn't run tests so a few things may be broken

 

Ondrej

--

You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

Mateusz

[Matthew Rocklin]

I like your use of the phrase "interpretation of Matrix Expressions." It makes it clear that we all probably have very different ideas for what Matrix Expressions should be. Hopefully there is a system which can handle a wide variety of needs. 

My current branch is here

https://github.com/mrocklin/sympy/tree/matrix_expr

The main difference between our two interpretations is that I think of MatrixExpr's as being built up mainly of MatrixSymbol's rather than Matrix objects. Probably there is a nice way to make both feasible. For my GSOC project I needed Block Matrices so those are better developed. There is a blogpost demonstrating them here.

[Aaron Meurer]

I think your implementations are completely different, though perhaps
they could be merged. Mateusz's implementation is more about allowing
unevaluated functions on Matrices, whereas Matthews is more about
expressions where the variables are to be understood as matrices (I
use this instead of the ambiguous "symbolic matrices").

Both have their advantages. For example, the laplace_expand() from
Mateusz's version doesn't make sense in Matthew's version. On the
other hand, Matthew's version truly lets you manipulate undefined
matrices, and can be a more efficient way of manipulating defined ones
(e.g., Identity(100000) vs. eye(100000)).

But other than some functions like laplace_expand() only making sense
for one, I don't see any reason why these two approaches can't or
shouldn't be merged. Actually, this is basically what was being
suggested when people said subs should work on symbolic matrices (a la
Matthew).

Aaron Meurer

[Mateusz Paprocki]

Hi,

On 30 July 2011 07:27, Aaron Meurer <asme...@gmail.com> wrote:

I think your implementations are completely different, though perhaps
they could be merged.  Mateusz's implementation is more about allowing
unevaluated functions on Matrices, whereas Matthews is more about
expressions where the variables are to be understood as matrices (I
use this instead of the ambiguous "symbolic matrices").

Both have their advantages.  For example, the laplace_expand() from
Mateusz's version doesn't make sense in Matthew's version.  On the
other hand, Matthew's version truly lets you manipulate undefined
matrices, and can be a more efficient way of manipulating defined ones
(e.g., Identity(100000) vs. eye(100000)).

But other than some functions like laplace_expand() only making sense
for one, I don't see any reason why these two approaches can't or
shouldn't be merged.  Actually, this is basically what was being
suggested when people said subs should work on symbolic matrices (a la
Matthew).

I'm sure that, sooner or later, those approaches will have to be merged, because those are really two views of a very similar (if not the same) problem domain. My original motivation came from reading lecture notes for undergraduates about the finite element method. As usually there was an introduction to basics of algebra needed to understand the later material, and my question was why it must be so hard to do it in SymPy (if possible at all). My branch is about "symbolic matrices" with explicit content. However, I don't see any problem with allowing transition between those two views (well at least in one direction). Suppose we have expr = Eq(A*x, b), where A, x, b are matrices/vectors of appropriate shape. First, I would like to be able to manipulate the expression alone, check various shapes (and ask SymPy if it makes sense), etc. Then I would like to write something like expr.expand(fullform=True) and get the same but with MatrixExpr with explicit indexed symbols or values (if entities like zeros or ones matrix was used in expr). Then I would like to make further transformation on this "full form".

Mateusz

[Aaron Meurer]

Also, Mateusz's MatrixExpr should be renamed to ImmutableMatrix.

Aaron Meurer

[Tim Lahey]

On Sat, Jul 30, 2011 at 2:31 AM, Mateusz Paprocki <mat...@gmail.com> wrote:

> I'm sure that, sooner or later, those approaches will have to be merged,
> because those are really two views of a very similar (if not the same)
> problem domain. My original motivation came from reading lecture notes
> for undergraduates about the finite element method. As usually there was an
> introduction to basics of algebra needed to understand the later material,
> and my question was why it must be so hard to do it in SymPy (if possible at
> all). My branch is about "symbolic matrices" with explicit content. However,
> I don't see any problem with allowing transition between those two views
> (well at least in one direction). Suppose we have expr = Eq(A*x, b), where
> A, x, b are matrices/vectors of appropriate shape. First, I would like to be
> able to manipulate the expression alone, check various shapes (and ask SymPy
> if it makes sense), etc. Then I would like to write something like
> expr.expand(fullform=True) and get the same but with MatrixExpr with
> explicit indexed symbols or values (if entities like zeros or ones matrix
> was used in expr). Then I would like to make further transformation on this
> "full form".

I would like to do things like differentiate c^T*A*c with respect to
the vector c. It's a common thing for finite elements. But I'd also
like block matrices. However, if you have symbolic matrix expressions,
you can put them in a matrix and then perform standard matrix
calculations on that and you'd have block matrix support. So, as long
as that's possible, there's no problem.

I've got by handling differentiating matrix-vector expressions in
Maple using their non-commutative support along with hacking together
handling the transpose and differentiation of it. But, I'd like proper
support.

If Sympy could support block matrices, that would be extremely useful
in control theory (where they're used all the time).

Cheers,

Tim.

--
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
http://about.me/tjlahey

[Andy Ray Terrel]

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
>
>

It might be a good idea to collect these uses in a wiki page that
descirbes the "User Stories", borrowing from the SCRUM process.

@Mateusz, doing symbolic FEM is not crazy there are lots of people who
have done it, including yours truly. The issue is that FEM is usually
done for very large problems for which symbolics is prohibatively
expensive. The FEniCS, Sundance, and others approach is to use
symbolics to describe operators for building a numeric matrix, SyFi
took your approach and build symbolics down to the matrix level but
for difficult problems had real problems with the expansion of the
symbolic expressions, Dohp is doing the symbolics at the quadrature
points in a mesh to evaluate the residuals, but again is having
trouble with the symbolic expressions.

-- Andy

[Matthew Rocklin]

I see that Andy just posted while I was writing. I'll post anyway although Ii thnk maybe the wiki page is a better start. 

It seems like we have a few people who want to contribute to the same concept. Should we put something on the master branch so that people can start adding to it? 

What would a general framework for matrix expressions look like that could handle most of the use-cases? What are the use cases? So far we have

• Expression that may contain matrices as atomic symbols
• Expression that may contain immutable matrices
• Derivatves of symbols
• Derivatves over indices?
• Block matrices

On Sat, Jul 30, 2011 at 1:40 AM, Tim Lahey <tim....@gmail.com> wrote:

--

[Matthew Rocklin]

Start of a wiki-page is here if people want to go this route. I put down the things that I think about. 

https://github.com/sympy/sympy/wiki/Matrix-Expressions

I'm also happy to continue the conversation over e-mail. 

What would be a good next step? How can I stimulate activity on this topic?

[Alan Bromborsky]

I looked at your wiki page  and since you teach linear algebra you might find the following book interesting -

http://www.amazon.com/Linear-Geometric-Algebra-Alan-Macdonald/dp/1453854932

also if you have not already done so you should look at the documentation for the geometric algebra module for sympy.

[Matthew Rocklin]

Will do!

As a disclaimer I don't actually teach Linear Algebra (though I respect those who do). The narratives posted are fictitious. I'm hypothesizing why people might want Matrix Expressions. 

[Aaron Meurer]

On Tue, Aug 2, 2011 at 9:00 AM, Matthew Rocklin <mroc...@gmail.com> wrote:
> Will do!
> As a disclaimer I don't actually teach Linear Algebra (though I respect
> those who do). The narratives posted are fictitious. I'm hypothesizing why
> people might want Matrix Expressions.

That's a little misleading. Maybe you should put them in the second person.

Aaron Meurer