0

I have a function that creates a 2-dim array, a Vandermonde matrix and is called as:

vandermonde(generator, rank)

Where generator is a n-sized array for example

generator = np.array([-1/2, 1/2, 3/2, 5/2, 7/2, 9/2])

and rank=4

Then I need to create 4 Vandermonde matrices (because rank=4) skewed by h in my space (that h is arbitrary here, lets call h=1).

Therefore I came with the following deterministic code:

V = np.array([
        vandermonde(generator-0*h, rank),
        vandermonde(generator-1*h, rank),
        vandermonde(generator-2*h, rank),
        vandermonde(generator-3*h, rank)
        ])

Then I want instead do multiple manual calls to vandermonde I used a for-loop as in:

V=[]
for i in range(rank):
    V.append(vandermonde(generator - h*i, rank))
V = np.array(V)

This approach works fine, but seems too "patchy". I tried a np.append approach as below:

M = np.array([])
for i in range(rank):
    M = np.append(M,[vandermonde(generator - h*i, rank)])

But didn't worked as I expected, seems np.append expand the array instead to create a new element.

My questions are:

  1. How can I not use standard Python lists, use directly a np approach cause np.append seems not behave as I expect, instead it just grow that array instead add a new array element

  2. Is there any more direct numpy approaches to this?

My vandermonde function is:

def vandermonde(generator, rank=None):
    """Returns a vandermonde matrix

    If rank not passwd returns a square vandermonde matrix
    """

    if rank is None:
        rank = len(generator)

    return np.tile(generator,(rank,1)) ** np.array(range(rank)).reshape((rank,1))

The expected answer is a 3 dimensional array with size (generator, rank, rank) where each element is one of the generator skewed vandermonde matrices. For the constants above(generator, rank, h) we have:

V= array([[[  1.  ,   1.  ,   1.  ,   1.  ,   1.  ,   1.  ],
        [ -0.5 ,   0.5 ,   1.5 ,   2.5 ,   3.5 ,   4.5 ],
        [  0.25,   0.25,   2.25,   6.25,  12.25,  20.25],
        [ -0.12,   0.12,   3.38,  15.62,  42.88,  91.12]],

       [[  1.  ,   1.  ,   1.  ,   1.  ,   1.  ,   1.  ],
        [ -1.5 ,  -0.5 ,   0.5 ,   1.5 ,   2.5 ,   3.5 ],
        [  2.25,   0.25,   0.25,   2.25,   6.25,  12.25],
        [ -3.38,  -0.12,   0.12,   3.38,  15.62,  42.88]],

       [[  1.  ,   1.  ,   1.  ,   1.  ,   1.  ,   1.  ],
        [ -2.5 ,  -1.5 ,  -0.5 ,   0.5 ,   1.5 ,   2.5 ],
        [  6.25,   2.25,   0.25,   0.25,   2.25,   6.25],
        [-15.62,  -3.38,  -0.12,   0.12,   3.38,  15.62]],

       [[  1.  ,   1.  ,   1.  ,   1.  ,   1.  ,   1.  ],
        [ -3.5 ,  -2.5 ,  -1.5 ,  -0.5 ,   0.5 ,   1.5 ],
        [ 12.25,   6.25,   2.25,   0.25,   0.25,   2.25],
        [-42.88, -15.62,  -3.38,  -0.12,   0.12,   3.38]]])

Some related ideas can be found in this discussion on: efficient-way-to-compute-the-vandermonde-matrix

Cœur
  • 37,241
  • 25
  • 195
  • 267
Lin
  • 1,145
  • 11
  • 28
  • Share implementation of `vandermonde`? – Divakar Jun 20 '18 at 12:51
  • @Divakar sure. Added. – Lin Jun 20 '18 at 12:53
  • ppl, just downvote for free? no comments added. – Lin Jun 20 '18 at 12:55
  • 1
    Numpy is efficient because it works with fix size arrays. The method `np.append` is no really efficient and is indeed expending an array (which means modification of the size => inneficient). To be efficient, you need to create an array at the right size before placing element into it. Since I do not completely get what the desired input / outputs are, it's hard to come up with a better proposition. – Mathieu Jun 20 '18 at 12:56
  • adding expected output and array dimensions. – Lin Jun 20 '18 at 12:58

1 Answers1

1

Use broadcasting to get the final 3D array in a vectorized manner -

r = np.arange(rank)
V_out = (generator - h*r[:,None,None]) ** r[:,None]

We can also use cumprod to achieve the exponential values for another solution -

gr = np.repeat(generator - h*r[:,None,None], rank, axis=1)
gr[:,0] = 1
out = gr.cumprod(1)
Divakar
  • 218,885
  • 19
  • 262
  • 358
  • Seems master broadcasts is a must have. I think that will even help me to rewrite the ```vandermonde``` function. – Lin Jun 20 '18 at 13:11
  • @Lin Yes, it will be : `(generator - h*i) ** np.arange(rank)[:,None]`. – Divakar Jun 20 '18 at 13:12
  • @Lin This could be helpful as well: https://stackoverflow.com/questions/48245987/efficient-way-to-compute-the-vandermonde-matrix – kmario23 Jun 20 '18 at 14:14
  • 1
    @kmario23 thank you, that really add some ideas to the problem. – Lin Jun 20 '18 at 14:51