views:

71

answers:

1

Hi, I have the following array in NumPy:

A = array([1, 2, 3])

How can I obtain the following matrices (without an explicit loop)?

B = [ 1 1 1
      2 2 2
      3 3 3 ]

C = [ 1 2 3
      1 2 3
      1 2 3 ]

Thanks!

+7  A: 

Edit2: The OP asks in the comments how to compute

n(i, j) = l(i, i) + l(j, j) - 2 * l(i, j)

I can think of two ways. I like this way because it generalizes easily:

import numpy as np

l=np.arange(9).reshape(3,3)
print(l)
# [[0 1 2]
#  [3 4 5]
#  [6 7 8]]

The idea is to use np.ogrid. This defines a list of two numpy arrays, one of shape (3,1) and one of shape (1,3):

grid=np.ogrid[0:3,0:3]
print(grid)
# [array([[0],
#        [1],
#        [2]]), array([[0, 1, 2]])]

grid[0] can be used as a proxy for the index i, and grid[1] can be used as a proxy for the index j.

So everywhere in the expression l(i, i) + l(j, j) - 2 * l(i, j), you simply replace i-->grid[0], and j-->grid[1], and numpy broadcasting takes care of the rest:

n=l[grid[0],grid[0]] + l[grid[1],grid[1]] + 2*l
print(n)
# [[ 0  6 12]
#  [10 16 22]
#  [20 26 32]]

However, in this particular case, since l(i,i) and l(j,j) are just the diagonal elements of l, you could do this instead:

d=np.diag(l)
print(d)
# [0 4 8]

d[np.newaxis,:] pumps up the shape of d to (1,3), and d[:,np.newaxis] pumps up the shape of d to (3,1).

Numpy broadcasting pumps up d[np.newaxis,:] and d[:,np.newaxis] to shape (3,3), copying values as appropriate.

n=d[np.newaxis,:] + d[:,np.newaxis] + 2*l
print(n)
# [[ 0  6 12]
#  [10 16 22]
#  [20 26 32]]

Edit1: Usually you do not need to form B or C. The purpose of Numpy broadcasting is to allow you to use A in place of B or C. If you show us how you plan to use B or C, we might be able to show you how to do the same with A and numpy broadcasting.


(Original answer):

In [11]: B=A.repeat(3).reshape(3,3)

In [12]: B
Out[12]: 
array([[1, 1, 1],
       [2, 2, 2],
       [3, 3, 3]])

In [13]: C=B.T

In [14]: C
Out[14]: 
array([[1, 2, 3],
       [1, 2, 3],
       [1, 2, 3]])

or

In [25]: C=np.tile(A,(3,1))

In [26]: C
Out[26]: 
array([[1, 2, 3],
       [1, 2, 3],
       [1, 2, 3]])

In [27]: B=C.T

In [28]: B
Out[28]: 
array([[1, 1, 1],
       [2, 2, 2],
       [3, 3, 3]])

From the dirty tricks department:

In [57]: np.lib.stride_tricks.as_strided(A,shape=(3,3),strides=(4,0))
Out[57]: 
array([[1, 1, 1],
       [2, 2, 2],
       [3, 3, 3]])

In [58]: np.lib.stride_tricks.as_strided(A,shape=(3,3),strides=(0,4))
Out[58]: 
array([[1, 2, 3],
       [1, 2, 3],
       [1, 2, 3]])

But note that these are views of A, not copies (as were the solutions above). Changing B, alters A:

In [59]: B=np.lib.stride_tricks.as_strided(A,shape=(3,3),strides=(4,0))

In [60]: B[0,0]=100

In [61]: A
Out[61]: array([100,   2,   3])

unutbu
I want to compute `n(i, j) = l(i, i) + l(j, j) - 2 * l(i, j)` for all i, j and get the result in matrix N
Yassin
Thanks a lot for your complete answer!
Yassin
++ for i,j = ogrid[...], a general method. (Have you thought of collecting your excellent posts, e.g. under scipy.org/Cookbook ?)
Denis
@Denis: I was all excited about contributing to http://www.scipy.org/Cookbook, but I'm stymied. It's taking me 133 seconds just to load the front page. Although there is a lot of great information there, and I'm not really sure mine is up to snuff, the site as a practical matter is just unworkable for me.
unutbu
@Unutbu, sad but true; have to look for another site. But the goal is important. How about a SO question along the lines "Suggestions and examples for organizing programming wikis":do you have suggestions, good examples, good practicesfor organizing programming wikis ?Most wikis are loose collections of specific Q+As;how can one organize such into coherent chapters with overview, tutorial, then Q+As ?(Yes this is programming related -- good programmers learn together.)
Denis