<html>
<head>
<BASE HREF="http://www.numeric-quest.com/haskell/Tensor.html">
<title>
Ndimensional tensors
</title>
</head>
<body>
<ul>
<center>
<h1>
***
</h1>
<h1>
Ndimensional tensors
</h1>
<b>
<br>
Literate Haskell module <i>Tensor.lhs</i>
</b>
<p>
Jan Skibinski, <a href="http://www.numeric-quest.com/news/">
Numeric Quest Inc.</a>, Huntsville, Ontario, Canada
<p>
1999.10.08, last modified 1999.10.16
</center>
<p>
<hr>
<p>
<i>
This is a quick sketch of what might be a basis of a real
Tensor module. This module has quite a few limitations (listed below).
I'd like to get some feedback on what should be a better
way to design it properly. Nevertheless, this module works
and is able to tackle complex and mundane manipulations
in the very straightforward way.
<p>
There are few arbitrary decisions we have taken. For example,
we consider a scalar to be a tensor of rank 0. This forces us to
do conversions between true scalars and such tensors, but it also
saves us a lot of headache related to typing restrictions. This
is a typical price paid for (too much?) generalization.
<p>
To get rid of those awful sums appearing in multiplications
of tensors we do introduce Einstein's summation convention by the way of
text examples
Hopefully it is clear and be well appreciated for its economy
of notation, which is standard in the tensor calculus.
<p>
Datatype <code>Tensor</code> defined here is an instance
of class <code>Eq</code>, <code>Show</code> and <code>Num</code>.
That means that one can compare tensors for equality and perform
basic numerical calculations, such as addition, negation,
subtraction, multiplication, etc.
<code>(==), (/=), (+), (), (*)</code>. In addition, several
customized operations, such as <code> (<*>)</code>
and <code>(<<*>>)</code> are defined for
variety of inner products.
<p>
Limitations of this module:
<ul>
<p>
<li>
Tensor components are Doubles. Why not Fraction, Complex, etc?
For a moment we will leave this question aside, and
return to it some time later. But we consider it
the important question
such generalization in some of our other modules:
<a href="http://www.numeric-quest.com/haskell/Orthogonals.html">
Orthogonals</a> and
<a href="http://www.numeric-quest.com/haskell/fractions.html">
Fraction</a>.
<p>
<li>
We are well aware that the decision to represent tensors
as nested objects will have significant impact on access
(and update
arrays seem to be better suited for such tasks, where all
indices must be explicitely computed first, but the access
time is linear. In contrary, the hierarchical data structure
defined here require very little effort in index computing
but the access time depends on the depth of the data tree.
<p>
But speed has not been tested yet, so we really do not know
how inefficient this module is and all of the above is
just a pure speculation. Certain operations of this module
seem to be quite well matched with this treelike data structure,
and because of it this design decision might be not so bad
after all.
<p>
<li>
The shape of tensors defined here involves two parameters:
dimension and rank. Rank is associated with the
depth of the tensor tree and corresponds to a total number
of indices by which you can access the individual components.
No limits are imposed on ranks and there are binary operations
which involve tensors of different ranks.
Dimension is associated with the breadth of the tree and
correspond to a number of values each index can take.
Dimension is fixed via constant <code>dims</code>. At first it might
seem as a severe limitation, but in fact one should never
mix tensors with different dimensions. One usually works
either with threedimensional tensors (classical mechanics,
electrodynamics, elasticity, etc.) or the fourdimentional
tensors (relativity theory).
</ul>
<p>
</i>
<p>
<hr>
<p>
<b>
Tensor datatype
</b>
<p>
<pre>
> module Tensor where
> import Data.Array(inRange)
> infixl 9 #
> infixl 9 ##
> infixl 7 <*>
> infixl 7 <<*>>
</pre>
Indices will assume values from range (1,dims) (defined below).
<p>
Tensor can contain a scalar value or a list of tensors.
This recursively defines tensor of any rank in nD space.
<pre>
> data Tensor = S Double
> | T [Tensor]
</pre>
There is no way we could specify the length of the list
<code>[Tensor]</code> in the data declaration. Typing is not
concerned with shapes.
We could of course use more specific representation of
this data structure, such as:
<pre>
data Tensor = S Double | T Tensor Tensor Tensor
</pre>
but then we would severily limit ourselves to threedimensional
tensors.
<p>
Rank is either 0 (scalars), 1 (vectors), or higher: 2, 3, 4 ...
<pre>
> rank :: Tensor -> Int
> rank t = rank' 0 t where
> rank' n (S _) = n
> rank' n (T xs) = rank' (n+1) (head xs)
</pre>
Here we define our tensor dimension as constant for this
module. All binary operations on tensors require the
same dimensions, so it makes sense to treat dimensions
as constants. But ranks can be different.
<pre>
> dims :: Int
> dims = 3
</pre>
<p>
<hr>
<p>
<b>
Showing
</b>
<p>
Tensors are printed as recursive lists with a word "Tensor"
prepended
<pre>
> instance Show Tensor where
> showsPrec 0 (S a) = showString "Tensor " . showsPrec 0 a
> showsPrec n (S a) = showsPrec n a
> showsPrec 0 (T xs) = showString "Tensor " . showList' 0 xs
> showsPrec n (T xs) = showList' n xs
> showList' :: (Show t) => Int -> [t] -> String -> String
> showList' _ [] = showString "[]"
> showList' n (x:xs) = showChar '[' . showsPrec (n+1) x . showRem (n+1) xs
> where
> showRem _ [] = showChar ']'
> showRem o (y:ys) = showChar ',' . showsPrec o y . showRem o ys
</pre>
<p>
<hr>
<p>
<b>
Input
</b>
<p>
Although tensors are printed as structured list
it is easier to input data via flat lists.
But make sure that the length of the list is one
of: dims^0, dims^1, dims^2, dims^3, dims^4, etc.
<p>
This function is quite inefficient for ranks higher than 4.
Compare, for example, timings of:
<pre>
tensor [1..3^6]
tensor [1..3^3] * tensor [1..3^3]
</pre>
Although both expressions create tensors of the same rank 6,
but the execution of the latter is much faster. This is
because the function <code>tensor</code> spends much
of its effort on recursively restructuring the flat lists
into the listsoflistsoflists...
<pre>
> tensor :: [Double] -> Tensor
> tensor xs
> | size == 1 = S (head xs)
> | q /= 0 = error "Length is not a power of dims"
> | otherwise = T (tlist p xs)
> where
> (p,q) = rnk 1 (quotRem size dims)
> rnk m (1, v) = (m, v)
> rnk m (u, 0) = rnk (m+1) (quotRem u dims)
> rnk m (_, v) = (m, v)
> size = length xs
> group n ys = group' n ys [] where
> group' o zs as
> | length zs == 0 = reverse as
> | length zs < o = reverse (zs:as)
> | otherwise = group' o (drop o zs) ((take o zs):as)
>
> tlist :: Int -> [Double] -> [Tensor]
> tlist 1 zs = map S zs
> tlist rnl zs = tlist' (rnl1) (map S zs)
> where
> tlist' 0 fs = fs
> tlist' o fs = tlist' (o1) $ map T $ group dims fs
</pre>
<p>
<hr>
<p>
<b>
Extraction and conversion
</b>
<p>
Tensor components are also tensors and can be extracted
via (#) operator
<pre>
> ( # ) :: Tensor -> Int -> Tensor
> (S a1) # 1 = S a1
> (S _) # _ = error "out of range"
> (T xs) # i = xs!!(i1)
> ( ## ) :: Tensor -> [Int] -> Tensor
> a ## [] = a
> a ## (x:xs) = (a#x) ## xs
</pre>
Tensors of rank 0 can be converted to scalars; i.e.,
simple numbers of type Double.
<pre>
> scalar :: Tensor -> Double
> scalar (S a) = a
> scalar (T _) = error "rank not 0"
</pre>
Tensors of rank 1 can be converted to vectors; i.e.,
lists with "dims" components of type Double
<pre>
> vector :: Tensor -> [Double]
> vector (S _) = error "rank not 1"
> vector a@(T xs)
> | rank a /= 1 = error "rank not 1"
> | otherwise = map scalar xs
</pre>
<p>
<hr>
<p>
<b>
Useful tensors: epsilon and delta
</b>
<p>
Function "epsilon' i j k" emulates values of the pseudotensor Eijk.
It is valid only for threedimensional tensors.
It takes three indices i,j,k from the range (1,3)
and returns one of the three values:
0.0, 1.0, 1.0
<pre>
> epsilon' :: Int -> Int -> Int -> Double
> epsilon' i j k
> | dims /= 3 = error "not 3-dims"
> | outside (1,3) i j k = error "Not in range"
> | (i == j) || (i == k) || (j == k) = 0
> | otherwise = epsilon1 i j k
> where
> epsilon1 m n o
> | (m == 1) && (n == 2) && (o == 3) = 1
> | (m == 3) && (n == 2) && (o == 1) = 1
> | otherwise = epsilon1 n o m
> outside (p,q) a b c =
> (not $ inRange (p,q) a) ||
> (not $ inRange (p,q) b) ||
> (not $ inRange (p,q) c)
</pre>
Function "delta' i j" emulates Kronecker's delta:
<pre>
> delta' :: Int -> Int -> Double
> delta' i j
> | i == j = 1
> | otherwise = 0
</pre>
Delta' and epsilon' can be converted to tensors
<pre>
> delta, epsilon :: Tensor
> delta = tensor [delta' i j | i <- [1..dims], j <- [1..dims]]
> epsilon = tensor [epsilon' i j k | i <- [1..3], j <- [1..3], k <- [1..3]]
</pre>
The components delta[ij] and epsilon[i,j,k] can be extracted
and converted to numbers. For example:
<pre>
scalar (epsilon#1#2#3) = 1
scalar (epsilon#1#1#3) = 0,
scalar (epsilon#3#2#1) = 1
</pre>
<p>
<hr>
<p>
<b>
Dot product
</b>
<p>
Dot product of two tensors of rank 1 could be defined as
tensor of rank 0. This is not the most efficient implementation
but we still want the dot product to be recognised as
tensor, so we loose on speed here:
<pre>
> dot :: Tensor -> Tensor -> Tensor
> dot a b = S (sum [scalar (a#i) * scalar (b#i) | i <- [1..dims]])
</pre>
<p>
<hr>
<p>
<b>
Cross product valid for 3D space only
</b>
<p>
The cross product of two vectors is another vector:
C = A x B. The pseudotensor Eijk is used to compute
such cross product.
<p>
First, here are numerical components of C, C[i]:
<pre>
> cross' :: Tensor -> Tensor -> Int -> Double
> cross' a b i = sum [(epsilon' i j k)* scalar (a#j) * scalar (b#k)|
> j<-[1..3],k<-[1..3], j/=k]
</pre>
And here is the full vector C (as tensor of rank 1):
<pre>
> cross :: Tensor -> Tensor -> Tensor
> cross a b = tensor (map (cross' a b) [1..3])
</pre>
Example:
<pre>
cross (tensor [1..3]) (tensor [1,8,1]) ==> Tensor [22.0, 2.0, 6.0]
</pre>
<p>
<hr>
<p>
<b>
Equality of tensors
</b>
<p>
Tensor can be admitted to class <code>Eq</code>. We only need to
define either equality or nonequality operation. We've chosen
to define the former: two tensors are equal if they have the same
rank and equal components:
<pre>
> instance Eq Tensor where
> (==) a b
> | ranka /= rank b = False
> | ranka == 0 = scalar a == scalar b
> | otherwise = and [(a#i) == (b#i) | i <- [1..dims]]
> where
> ranka = rank a
>
</pre>
<p>
<hr>
<p>
<b>
Tensor as instance of class Num
</b>
<p>
To admit tensors to class <code>Num</code> we have to
support all the operations from that class. Here is
the class Num declaration taken from the Prelude:
<pre>
class (Eq a, Show a) => Num a where
(+), (), (*) :: a -> a -> a
negate :: a -> a
abs, signum :: a -> a
fromInteger :: Integer -> a
x y = x + negate y
negate x = 0 x
</pre>
All operations but <code>(*)</code> are straightforward,
meaningful and easy to implement. The semantics of multiplication
<code>(*)</code> is, however, not so obvious and it is up to us
how to define it: as an inner product or as an outer
product. We have chosen the latter, which means that the
operation <code>c = a * b</code> produces a new tensor <code>c</code>
whose rank is a sum of the ranks of tensors being
multiplied:
<pre>
rank c = rank a + rank b
</pre>
Suffice to add that tensor products are generally not
commutative; that is:
<pre>
a * b /= b * a
</pre>
That said, here is the instantiation of <code>Num</code>
for datatype Tensor:
<pre>
> instance Num Tensor where
> (+) a b
> | ranka /= rank b = error "different ranks"
> | ranka == 0 = S (scalar a + scalar b)
> | otherwise = T [a#i + b#i | i <- [1..dims]]
> where
> ranka = rank a
> negate (S a1) = S (negate a1)
> negate (T xs) = T (map negate xs)
> abs (S a1) = S (abs a1)
> abs (T xs) = T (map abs xs)
> signum (S a1) = S (signum a1)
> signum (T xs) = T (map signum xs)
> fromInteger n = S (fromInteger n)
> (*) (S a1) (S b1) = S (a1*b1)
> (*) a@(S _) (T xs) = T (map (a*) (take dims xs))
> (*) (T xs) b = T (map (*b) (take dims xs))
</pre>
Having defined the operation <code>(*)</code> as an outer product
such operation will generally increase the rank of the outcome.
For example, if <code>a</code> is a tensor of rank 2 (matrix) and
<code>b</code> is a tensor of rank 1 (vector) then the result is
a tensor of rank 3:
<pre>
c = a * b, that is
c[ijk] = a[ij] b[k]
</pre>
But this is not what is typically considered a multiplication
of tensors; we are more often than not interested in the inner
products, informally described below.
<p>
<hr>
<p>
<b>
Contraction
</b>
<p>
<p>
Eistein's indexing convention of tensors is based on
the distinction between free indices and bound indices.
Free indices appear in the tensorial expressions, such
as <code>A[ijkl]</code>, once only and they indicate
a freedom for substitution of any specific index
from the range of valid indices. This range is (1,3)
for 3D tensors. The expression <code>A[ijkl]</code>
represents in fact one of 3^4 possible components
of the tensor <code>A</code