0

Consider a NumPy ndarray A of floats, with a dimension of n and an arbitrary shape of D=[d1,...,dn] (dis are nonnegative integers). How can I populate A to have for example:

A[j1,...,jn]=sqrt(j1*...*jn)

where 0<=ji<=di. If I knew n and it was fixed then I could nest n for loops to simply populate the ndarray. however in my program that's not the case, apart from the fact that will not be efficient. I was wondering if there is any way to

  • populate the ndarray A of the given arbitrary shape D (for example with the formula above or any other nonthrivial function of the indices).
  • preferably avoid using python for loops and take advantage of NumPy's underlying functionality.

Thanks for your help in advance.

Foad S. Farimani
  • 12,396
  • 15
  • 78
  • 193
  • Do I understand correctly that you want to populate your randomly-shaped array with (in this example) the square root of the product of the indices? In particular, that *the product of the indices* is the part you're struggling with? – 9769953 Aug 10 '18 at 08:34
  • @9769953 mostly true. The shape and dimension are arbitrary given by `D`. the function of indices is also not important. it can be anything, not necessarily the product. my final goal is to implement what I have suggested [here](https://math.stackexchange.com/questions/2877478/cauchy-product-of-multivariate-formal-power-series). – Foad S. Farimani Aug 10 '18 at 08:38

2 Answers2

0

One important fact to realize is that you can use broadcasting to solve this problem efficiently. So for the 2D case you could do

d1, d2 = (3, 4)
A = numpy.sqrt(numpy.arange(d1)[:, None] * numpy.arange(d2)[None, :])
# array([[0.        , 0.        , 0.        , 0.        ],
#        [0.        , 1.        , 1.41421356, 1.73205081],
#        [0.        , 1.41421356, 2.        , 2.44948974]])

Once you feel comfortable with using broadcasting to do these outer products (or sums, or comparisons etc.) we can try to solve this for the nD case.

Looking at the input arrays of the above code we realize they have shapes

(d1,  1)
( 1, d2)

So to do this in nD we need to find a method that takes linear index arrays and automatically creates arrays of shapes

(d1,  1,  1, ...)
( 1, d2,  1, ...)
( 1,  1, d3, ...)

Numpy offers such a function: numpy.meshgrid(..., sparse=True)

numpy.meshgrid(numpy.arange(3), numpy.arange(4), sparse=True)

Knowing this we can put it all together in one line:

D = (3, 4, 5)
numpy.sqrt(numpy.prod(numpy.meshgrid(*[numpy.arange(d) for d in D], sparse=True, indexing='ij')))
# array([[[0.        , 0.        , 0.        , 0.        , 0.        ],
#         [0.        , 0.        , 0.        , 0.        , 0.        ],
#         [0.        , 0.        , 0.        , 0.        , 0.        ]],
# 
#        [[0.        , 0.        , 0.        , 0.        , 0.        ],
#         [0.        , 1.        , 1.41421356, 1.73205081, 2.        ],
#         [0.        , 1.41421356, 2.        , 2.44948974, 2.82842712]],
# 
#        [[0.        , 0.        , 0.        , 0.        , 0.        ],
#         [0.        , 1.41421356, 2.        , 2.44948974, 2.82842712],
#         [0.        , 2.        , 2.82842712, 3.46410162, 4.        ]],
# 
#        [[0.        , 0.        , 0.        , 0.        , 0.        ],
#         [0.        , 1.73205081, 2.44948974, 3.        , 3.46410162],
#         [0.        , 2.44948974, 3.46410162, 4.24264069, 4.89897949]]])

Performance evaluation

To evaluate the performance of all three solutions, let's time their speed for several different tensor sizes:

D=(2,3,4,5)

%timeit np.fromfunction(function=myfunc2, shape=D)
# 501 µs ± 9.34 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit np.fromfunction(function=creation_function, shape=D)
# 24.2 µs ± 455 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit numpy.sqrt(numpy.prod(numpy.meshgrid(*[numpy.arange(d) for d in D], sparse=True, indexing='ij')))
# 30.9 µs ± 1.02 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

D=(20,30,40,50)

%timeit np.fromfunction(function=myfunc2, shape=D)
# 4.64 s ± 36.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit np.fromfunction(function=creation_function, shape=D)
# 36.7 ms ± 1.17 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit numpy.sqrt(numpy.prod(numpy.meshgrid(*[numpy.arange(d) for d in D], sparse=True, indexing='ij')))
# 9 ms ± 237 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

D=(200,30,40,50)

%timeit np.fromfunction(function=myfunc2, shape=D)
# never completed
%timeit np.fromfunction(function=creation_function, shape=D)
# 508 ms ± 7.41 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit numpy.sqrt(numpy.prod(numpy.meshgrid(*[numpy.arange(d) for d in D], sparse=True, indexing='ij')))
# 88.1 ms ± 1.63 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

D=(200,300,40,50)

%timeit np.fromfunction(function=myfunc2, shape=D)
# never completed
%timeit np.fromfunction(function=creation_function, shape=D)
# 5.8 s ± 565 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit numpy.sqrt(numpy.prod(numpy.meshgrid(*[numpy.arange(d) for d in D], sparse=True, indexing='ij')))
# 1.29 s ± 15.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Nils Werner
  • 34,832
  • 7
  • 76
  • 98
  • Thanks a lot. This solution is great for the specific example provided in the OP. However, I am not sure if it is possible to extend it to an arbitrary function of indices as mentioned in the first requirement above. My final goal is to implement what I have suggested [here](https://math.stackexchange.com/questions/2877478/cauchy-product-of-multivariate-formal-power-series) – Foad S. Farimani Aug 10 '18 at 09:04
  • 1
    Well, if this answers the original question, mark it as answered and open another question asking how to implement that specific function (preferrably also providing some `for` loop based code, not everybody is able to quickly understand those formulas) – Nils Werner Aug 10 '18 at 09:10
  • I'm reading the `numpy.meshgrid` I will try to modify your code and see if I should ask a new question or accept any of the answers here. – Foad S. Farimani Aug 10 '18 at 09:13
  • That's fine, but don't move the goalpost. You asked *this question*, and *this answers it*. It's no problem if you have another question, just post another question. But don't go back and edit your original post to ask something different (rendering all of the answers invalid and confusing future readers). – Nils Werner Aug 10 '18 at 09:16
  • Try larger shapes or higher dimensions, it will eventually be faster. – Nils Werner Aug 10 '18 at 09:20
  • @NilsWerner thanks for your help. The other answer was acceptable because it satisfied the arbitrary function requirements. I appreciate your support. upvote. – Foad S. Farimani Aug 10 '18 at 09:22
  • It definitely allows arbitrary functions (within the limits of what vectorization allows), you just need to specify what function you want. – Nils Werner Aug 10 '18 at 09:23
  • as mentioned the other solution doesn't work for dimensions above 3! so I had to remove the acceptance flag. – Foad S. Farimani Aug 10 '18 at 09:28
  • I'm going try modifying your code in a way to take an arbitrary function of the indiciies. BTW thanks a lot for the comparison. your solution is definitely the fastest. – Foad S. Farimani Aug 10 '18 at 10:44
  • would you please compare your result against `np.fromfunction(lambda i, j, k, l: np.sqrt(i*j*k*l), D)` it seems to yield a slightly different result! – Foad S. Farimani Aug 10 '18 at 10:51
  • 1
    The two first dimensions were swapped. Setting the parameter `indexing='ij'` for `meshgrid()` fixed it, they're equal now. – Nils Werner Aug 10 '18 at 10:57
  • I think the reason your code is way faster is because of the `sparse=True` option. it takes advantage of the fact that in this specific case many elements of the ndarray are zero. which is not necessarily true for any arbitrary function of indices. would you please try the `sparse=False` option and see if the difference is still significant? – Foad S. Farimani Aug 10 '18 at 11:02
  • 1
    Your assumption is incorrect. It just creates an [open grid](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.ogrid.html), instead of a [full grid with millions of duplicate values](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.mgrid.html). – Nils Werner Aug 10 '18 at 11:09
  • I think I know have a rough understanding of what `numpy.meshgrid` with `sparse=True`. Now I can create a python function and then vectorize it to accept the output of `numpy.meshgrid`, or use a lambda function to do so. – Foad S. Farimani Aug 10 '18 at 11:30
  • 1
    Yup. Also, note that I am really intrigued to see how the problem from the other question translates into a `for` loop solution, and how one could vectorize it. Can you please create such a question now? :-D – Nils Werner Aug 10 '18 at 11:34
  • :)) I am starting the implementation now. I have to implement some functions: 1. slicing an nd array given a D, if any of the elements of D is bigger than the arrays dimension the rest will be zero 2. flipping the ndarray in all possible directions 3. inner product of two ndarrays ... – Foad S. Farimani Aug 10 '18 at 11:38
  • as promised, [here](https://gist.github.com/Foadsf/80abba07ea892e50da3df7e318844f33) the function. Now I need to vectorize this and populate the final array – Foad S. Farimani Aug 10 '18 at 19:38
  • I also posted a question [here](https://stackoverflow.com/questions/51794274/convolution-of-numpy-arrays-of-arbitrary-dimension-for-cauchy-product-of-multiva) as my attempts weren't successful so far. – Foad S. Farimani Aug 10 '18 at 21:32
  • Did this answer help you to solve your problem? If so, please remember to [mark it as accepted](https://meta.stackexchange.com/questions/5234/how-does-accepting-an-answer-work). – Nils Werner Aug 20 '18 at 09:00
  • It indeed helped and I upvoted it for that reason. But still does not satisfy all the criteria mentioned in the OP. The answer I have written below does satisfy but it is very inefficient. If you don't mind I want to leave it open. – Foad S. Farimani Aug 20 '18 at 09:03
  • It absolutely and 100% fulfils what was asked in the question. You *really* need to learn how to ask concise questions. – Nils Werner Aug 20 '18 at 09:10
  • you are absolutely write. I need to learn a lot regarding concise question writing. as discussed above, you have indeed implemented a highly efficient version of the example provided in the OP. However I had mentioned `nontrivial function of the indices`. My final goal is to use conv function mentioned [here ](https://stackoverflow.com/questions/51821540/how-to-implement-a-multidimensional-version-of-nympy-convolve). to populate the final ndarray. – Foad S. Farimani Aug 20 '18 at 09:15
  • But that nontrivial function **is not part of this question**, which makes this answer a correct one. – Nils Werner Aug 20 '18 at 09:16
  • Dear Nils, I really appreciate your support. Not accepting your result don't mean otherwise. the term `nontrivial` was mentioned in the OP. I accepted your post to show my gratitude. But I hope I succeed to convey the message. – Foad S. Farimani Aug 20 '18 at 09:20
-1

A modified version of ibarronds solution without using labmda and works with higher dimensions above 3:

import numpy as np

def myfunc(*J):
    return np.sqrt(np.prod(np.array(J)))
myfunc2=np.vectorize(myfunc)

D=(2,3,4,5)

np.fromfunction(function=myfunc2 , shape=D)

P.S. Unfortunately he has removed his answer so I copy it here for reference:

creation_function = lambda *args: np.sqrt(np.prod(np.array([*args]), axis=0))
np.fromfunction(function=creation_function, shape=D)
Foad S. Farimani
  • 12,396
  • 15
  • 78
  • 193
  • Good to hear that it worked out in the end! Mine worked too with a small tweak. Oh and lambda functions can also be turned into any function in my solution! – ibarrond Aug 10 '18 at 10:21
  • @ibarrond thanks a lot. your solution helped a lot. but I'm afraid it doesn't work. I left a comment. – Foad S. Farimani Aug 10 '18 at 10:23
  • Be aware that your solution is several orders of magnitude slower than the other two solutions. – Nils Werner Aug 10 '18 at 10:23
  • @NilsWerner the issue is ibarrond solution still doesn't work for me and yours is too specific for the example above. My final goal is to write a function which depends on the indices and then use that to populate the ndarray. yours first creates an ndarray of the given size and then applies numpy vectorized functions on the existing elements. – Foad S. Farimani Aug 10 '18 at 10:26
  • Why don't you just specify *what function you want to actually compute*, then? – Nils Werner Aug 10 '18 at 10:27
  • Anyhow, @Foad you ended up solving your own question. Set it as right for future references! – ibarrond Aug 10 '18 at 10:28
  • @NilsWerner the function at the end of [this question](https://math.stackexchange.com/questions/2877478/cauchy-product-of-multivariate-formal-power-series) is what I'm trying to implement. your solution works but does not satisfy the requirement of `nontrivial function of the indices` as mentioned in the OP. – Foad S. Farimani Aug 10 '18 at 10:28
  • Sorry, but I don't understand those formulas. Can you translate it to some nested `for` loop code? – Nils Werner Aug 10 '18 at 10:30
  • @NilsWerner If your question is regarding my formula [here](https://math.stackexchange.com/questions/2877478/cauchy-product-of-multivariate-formal-power-series), then I'm working on an implementation now. – Foad S. Farimani Aug 10 '18 at 10:32
  • @ibarrond your code doesn't work with [this online interpreter](https://www.tutorialspoint.com/execute_python_online.php). – Foad S. Farimani Aug 10 '18 at 10:33