Subject: [UFL-dev] The role of * in UFL

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 7 Apr 2008 11:51:24 +0200

I'd like to hear opinions of how * should behave in different contexts.

To begin with, we will of course agree on the meaning of * between
scalar expressions, which may or may not have free indices. Repeated
free indices will be interpreted as implicit summation. Single free
indices are combined to define a larger set of free indices. This is
already implemented. (Some of the implicit summation is not, but I
have it on paper).

We may perhaps also agree that multiplying a matrix (rank 2 tensor)
with a vector (rank 1 tensor) means contraction over the last index of
the matrix and the index of the vector, which is equivalent with the
dot product. This will allow f.ex. Poissons equation with a conduction
tensor (M Du : Dv) to be written like inner(M*grad(u), grad(v)). Ok?

But if so, is this to be a special case, or do we simply define * to be dot?
This will eventually lead to (let u,v be vectors, A,B,C matrices, a scalar):

Very natural notation:
u = A*v = dot(A, v) = "A v" in linear algebra notation
C = A*B = dot(A, B) = "A B" in linear algebra notation

Feels natural to me:
u = v*A = dot(v, A) = "v^T A" in linear algebra notation
a = u*A*v = (u*A)*v = u*(A*v)
  = dot(u, A*v) = dot(u*A, v)
  = dot(dot(u, A), v) = dot(u, dot(A, v)) = "u^T A v" in linear algebra notation

However, keeping * and dot separate may make things simpler in other
ways in the implementation. And easy code understanding isn't always
the same as shorter expressions. So what do you think?

Note that dot(f, g) is defined differently in UFL than in FFC:
    dot(f, g) == sum(f[...,k] * g[k, ...] for k in range(dim))
where ... means the obvious for a tensor of rank > 1.
Using ... to "skip indices" is also implemented.
The contraction A:B is instead implemented as inner(A, B).

Another example, Navier-Stokes with explicit convection term (freely
from memory):
a = u - w + dt * ( w * grad(u) * v - nu*inner(grad(u), grad(v)) +
p*div(v) - f*v )
vs.
a = u - w + dt * ( dot(dot(w, grad(u)), v) - nu*inner(grad(u),
grad(v)) + p*div(v) - dot(f, v) )

--
Martin

---

From: Kent-Andre Mardal <[email-masked]>
Date: Mon, 07 Apr 2008 12:28:51 +0200

Difficult subject, in particular when it comes to matrices. 
I think we should let * and dot be the same, I mean dot is just a word
for *, right ? 

The meaning of *should be the same as in Matlab, i.e. 
A*B is not the contraction, but the matrix-matrix product. 

Kent

man, 07.04.2008 kl. 11.51 +0200, skrev Martin Sandve Alnæs:
> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 7 Apr 2008 12:57:14 +0200

>  Difficult subject, in particular when it comes to matrices.

I'd say it's most difficult when it comes to _vectors_, not matrices.
Actually, we're not dealing with vectors and matrices, we're dealing
with tensors. A rank 1 tensor is not the same as a vector in linear
algebra, there is no concept of rows or columns.

>  I think we should let * and dot be the same, I mean dot is just a word

Not really, that's the whole point. Did you look at the examples? "*"
means multiplication, but what does multiplication mean for tensors of
rank > 0? It can be interpreted differently, and that's why I'm asking
here.

>  A*B is not the contraction, but the matrix-matrix product.

That's one of the trivial examples in my previous email. A*B and A*v
work out fine, but v*A and u*v would be invalid products in
Matlab/linear algebra notation. If you write (with tensors) "u v", it
(often, usually, always?) means outer(u, v). Translating from linear
algebra notation to UFL tensor notation, "u v" is invalid, "u^T v" is
dot(u,v) and "u v^T" is outer(u, v).

"Paper notation" contains a lot of special cases that doesn't
generalize, and would lead to a tangled mess of code. That's why we
must make our own precise definitions, that are always well defined
and translates well to code. And since we're dealing with tensor
algebra, we must use tensor algebra as our basis and not linear
algebra where they differ.

-- 
Martin

---

From: Pearu Peterson <[email-masked]>
Date: Mon, 07 Apr 2008 13:12:01 +0200

> algebra where they differ.

FYI, in sympycore we have the same problem with matrices
and there we have resolved it as follows (for details see 
http://code.google.com/p/sympycore/wiki/MatrixSupportIdeas):

0) currently vectors are represented as 1-column matrices
1) by default, `A * B` is matrix product.
2) matrices have attributes .A, .T, etc (`.T` would be `^T` in your 
case) that return "array", "transpose", etc *views* of matrices
and these views define what kind of multiplication will be used.

For example, `A.A * B` would be equivalent to elementwise multiplication
of matrices. `A.T * B` would be dot product of matrices.
Outer products could be modelled as `A.O * B`, for instance.

Regards,
Pearu

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 7 Apr 2008 13:49:08 +0200

>  Outer products could be modelled as `A.O * B`, for instance.

I've already implemented A.T the same way. But I think f.ex. outer(A,
B) is clearer than A.O*B, and "A.O" in itself doesn't make sense to
me. Also, I've done inv(A) instead of A.I.

A.A*B we would currently have to do with manual loops, but I'm not
sure we need it.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 7 Apr 2008 14:26:31 +0200

> A.A*B we would currently have to do with manual loops, but I'm not

My suggestion would be to first have some very generic operators like                                                        
dot() and inner() which work for all kinds of tensors (when the                                                              
dimensions match). Then we add the * operator for a few linear algebra                                                       
operations that we want to support:                                                                                          
                                                                                                                             
  A*B      rank(A) == rank(B) == 2                                                                                           
  A*x      rank(A) == 2 and rank(x) == 1                                                                                     
  x.T*y    rank(x) == rank(y) == 1 and transposed(x) and not transposed(y)                                                   
  x.T*A    rank(x) == 1 and rank(A) == 2 and transposed(x)                                                                   
                                                                                                                             
Are there any others?                                                                                                        
                                                                                                                             
For outer products, one would need to write outer(x, y), which is very                                                       
natural.                                                                                                                     
                                                                                                                             
--                                                                                                                           
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 7 Apr 2008 14:35:09 +0200

>  For outer products, one would need to write outer(x, y), which is very

I don't agree with the x.T stuff. There is no differentiation between
columns and rows in tensor algebra, so the transpose has no meaning
for a rank 1 tensor. This has to generalize to where x is a rank 1
expression and not only an explicitly constructed vector. If you have
a rank 4 tensor and pick a sub-tensor of rank 1 from it, is it a
column or row? Columns and rows doesn't generalize, and special cases
are evil.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 7 Apr 2008 14:40:27 +0200

> column or row? Columns and rows doesn't generalize, and special cases

In that case, I don't think we should support x.T*y.

So we end up with just these operations:

  u*foo    rank(u) == 0
  foo*u    rank(u) == 0

In all other cases, one needs to use dot, inner or outer.

-- 
Anders

---

From: Pearu Peterson <[email-masked]>
Date: Mon, 07 Apr 2008 14:42:54 +0200

> column or row? Columns and rows doesn't generalize, and special cases

Transpose is well-defined also for arbitrary rank tensors. It swaps
the first two indices:

   transpose(t_ijkl..) -> t_jikl..

One can also defined generalised transpose:

   transpose(t_ijkl.., [p_i, p_j, p_k, p_l, ..]) -> t_(p_i)(p_j)(p_k)(p_l)..

Pearu

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 7 Apr 2008 14:54:07 +0200

> In all other cases, one needs to use dot, inner or outer.

I'm ok with that. In fact, that's the set of assertions I have
implemented in the current code.

In general, keeping the number of special definitions down reduces the
chance of misinterpretation.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 7 Apr 2008 14:56:43 +0200

> chance of misinterpretation.

ok!

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 7 Apr 2008 14:59:53 +0200

>    transpose(t_ijkl.., [p_i, p_j, p_k, p_l, ..]) -> t_(p_i)(p_j)(p_k)(p_l)..

We _could_ support something like that, but I don't think it's worth
the trouble since the more general index notation covers this without
making detailed definitions:

# reference components
Aijkl = A[i,j,k,l]
# make tensor with different axis ordering
B = Tensor(Aijkl, (j,i,k,l))

--
Martin

---

From: Pearu Peterson <[email-masked]>
Date: Mon, 07 Apr 2008 15:01:34 +0200

> B = Tensor(Aijkl, (j,i,k,l))

Nice:)

Pearu

---

From: Pearu Peterson <[email-masked]>
Date: Mon, 07 Apr 2008 16:03:04 +0200

The following thread might be of interest:
http://groups.google.com/group/sage-support/browse_thread/thread/cee2d7b5f889aa4b

Pearu

---

From: Anders Logg <[email-masked]>
Date: Mon, 7 Apr 2008 20:44:52 +0200

> http://groups.google.com/group/sage-support/browse_thread/thread/cee2d7b5f889aa4b

Interesting discussion. :-)

Some quotes:

  "It works nice, but it's quite difficult to install."

  "I tried Fenics about a year ago, so they may have improved."

-- 
Anders