[Numpy-discussion] Toward release 1.0 of NumPy
Charles R Harris
charlesr.harris at gmail.com
Thu Apr 13 17:14:32 EDT 2006
On 4/13/06, Tim Hochberg <tim.hochberg at cox.net> wrote:
>
> Charles R Harris wrote:
>
> >
> >
> > On 4/13/06, *Tim Hochberg* <tim.hochberg at cox.net
> > <mailto:tim.hochberg at cox.net>> wrote:
> >
> > Charles R Harris wrote:
> >
> > > Tim,
> > >
> > > On 4/13/06, *Tim Hochberg* <tim.hochberg at cox.net
> > <mailto:tim.hochberg at cox.net>
> > > <mailto:tim.hochberg at cox.net <mailto:tim.hochberg at cox.net>>>
> wrote:
> > >
> > > Alan G Isaac wrote:
> > >
> > > >On Thu, 13 Apr 2006, Charles R Harris apparently wrote:
> > > >
> > > >
> > > >>The Kronecker product (aka Tensor product) of two
> > > >>matrices isn't a matrix.
> > > >>
> > > >>
> > > >
> > > >That is an unusual way to describe things in
> > > >the world of econometrics. Here is a more
> > > >common way:
> > > >http://planetmath.org/encyclopedia/KroneckerProduct.html
> > > < http://planetmath.org/encyclopedia/KroneckerProduct.html>
> > > >I share Sven's expectation.
> > > >
> > > >
> > > mathworld also agrees with you. As does the documentation
> > (as best
> > > as I
> > > can tell) and the actual output of kron. I think Charles must
> be
> > > thinking of the tensor product instead.
> > >
> > >
> > > It *is* the tensor product, A \tensor B, but it is not the most
> > > general tensor with four indices just as a bivector is not the
> most
> > > general tensor with two indices. Numerically, kron chooses to
> > > represent the tensor product of two vector spaces a, b with
> > dimensions
> > > n,m respectively as the direct sum of n copies of b, and
> the tensor
> > > product of two operators takes the given form. More generally, the
> B
> > > matrix in each spot could be replaced with an arbitrary matrix
> > of the
> > > correct dimensions and you would recover the general tensor with
> > four
> > > indices.
> > >
> > > Anyway, it sounds like you are proposing that the tensor (outer)
> > > product of two matrices be reshaped to run over two indices. It
> > seems
> > > that likewise the tensor (outer) product of two vectors should be
> > > reshaped to run over one index ( i.e. flat). That would do the
> > trick.
> >
> > I'm not proposing anything. I don't care at all what kron does. I
> > just
> > want to fix the return type if that's feasible so that people stop
> > complaining about it. As far as I can tell, kron already returns a
> > flattened tensor product of some sort. I believe the general tensor
> > product that you are talking about is already covered by
> > multiply.outer,
> > but I'm not sure so correct me if I'm wrong. Here's what kron does
> as
> > present:
> >
> > >>> a
> > array([[1, 1],
> > [1, 1]])
> > >>> kron(a,a) # => 4x4 matrix
> > array([[1, 1, 1, 1],
> > [1, 1, 1, 1],
> > [1, 1, 1, 1],
> > [1, 1, 1, 1]])
> >
> >
> > Good at first look. Lets see a simpler version... Nevermind, seems
> > numpy isn't working on this machine (X86_64, fc5 64 bit) at the
> > moment, maybe I need to check out a clean version.
> >
> > >>> kron(a,a[0]) => 8x1
> > array([1, 1, 1, 1, 1, 1, 1, 1])
> >
> >
> > Looks broken. a[0] should be an operator (matrix), so either it should
> > be (2,1) or (1,2).
>
> Since a is an array here, a[0] is shape (2,). Let's repeat this
> excercise using matrices, which are always rank-2 and see if they make
> sense.
>
> >>> m
> matrix([[1, 1],
> [1, 1]])
> >>> kron(m, m[0])
> matrix([[1, 1, 1, 1],
> [1, 1, 1, 1]])
> >>> kron(m,m[:,0])
> matrix([[1, 1],
> [1, 1],
> [1, 1],
> [1, 1]])
>
> That looks OK.
>
> > In the first case, the return should have shape (4,2), in the latter
> > (2,4). Should probably raise an error as the result strikes me as
> > ambiguous. But I have to admit I am not sure what the point of this
> > particular construction is.
> >
> > >>> kron(a[0], a[0])
> > Traceback (most recent call last):
> > File "<stdin>", line 1, in ?
> > File "C:\Python24\Lib\site-packages\numpy\lib\shape_base.py", line
> > 577, in kron
> > result = concatenate(concatenate(o, axis=1), axis=1)
> > ValueError: 0-d arrays can't be concatenated
> >
> >
> >>> kron(m[0], m[0])
> matrix([[1, 1, 1, 1]])
> >>> kron(m[:,0], m[:,0])
> matrix([[1],
> [1],
> [1],
> [1]])
> >>> kron(m[:,0],m[0])
> matrix([[1, 1],
> [1, 1]])
>
> > See above. this could be (1,4) or (4,1), depending.
>
> All of these look like they're probably right without thinking about it
> too hard.
>
> >
> > >>> b.shape
> > (2, 2, 2)
> > >>> kron(b,b).shape
> > (4, 4, 2, 2)
> >
> >
> > I think this is doing transpose(outer(b,b), axis=(0,2,1,3)) and
> > reshaping the first 4 indices into 2. Again, I am not sure what the
> > point is for these operators. Now another way to get all this
> > functionality is to have a contraction function or method with a list
> > of axis. For instance, consider the matrices A(i,j) and B(k,l)
> > operating on x(j) and y(l) like A(i,j)x(j) and B(k,l)y(l), then the
> > outer product of all of these is
> >
> > A(i,j)B(k,l)x(j)y(l)
> >
> > with the summation convention on the indices j and l. The result
> > should be the same as kron(A,B)*kron(x,y) up to a permutation of rows
> > and columes. It is just a question of which basis is used and how the
> > elements are indexed.
> >
> > So, it looks like the 2d x 2d product obeys Alan's definition. The
> > other
> > products are probably all broken.
> >
> Here's my best guess as to what is going on:
> 1. There is a relatively large group of people who use Kronecker
> product as Alan does (probably the matrix as opposed to tensor math
> folks). I'm guessing it's a large group since they manage to write the
> definitions at both mathworld and planetmath.
> 2. kron was meant to implement this.
> 2.5 People who need the other meaning of kron can just use outer, so
> no real conflict.
> 3. The implementation was either inappropriately generalized or it
> was assumed that all inputs would be matrices (and hence rank-2).
Uh-huh.
Assuming 3. is correct, and I'd like to hear from people if they think
> that the behaviour in the non rank-2 cases is sensible, the next
> question is whether the behaviour in the rank-2 cases makes sense. It
> seem to, but I'm not a user of kron. If both of the preceeding are true,
> it seems like a complete fix entails the following two things:
> 1. Forbid arguments that are not rank-2. This allows all matrices,
> which is really the main target here I think.
> 2. Fix the return type issue. I have a fix for this ready to commit,
> but I want to figure out the first part as well.
I think it was inappropriately generalized, it is hard to make sense of what
kron means for rank > 2. So I vote for restricting the usage to matrices, or
arrays of rank two. This avoids the both the ambiguity of rank one arrays
and big why that arises for arrays with rank > 2. Note that in tensor
algebra the rank 1 problem is solved by the use of upper or lower indices,
lower index => [1,n], upper index => [n,1].
Hmm, I should to check that kron is associative: kron(kron(a,b),c) ==
kron(a, kron(b,c)) like a good tensor product should be. I suspect it is.
Regards,
>
> -tim
Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20060413/25f1192f/attachment-0001.html>
More information about the NumPy-Discussion
mailing list