0

Currently I have a C matrix generated by:

def c_matrix(n):
    exp = np.exp(1j*np.pi/n)
    exp_n = np.array([[exp, 0], [0, exp.conj()]], dtype=complex)
    c_matrix = np.array([exp_n**i for i in range(1, n, 1)], dtype=complex)
    return c_matrix

What this does is basically generate a list of number from 0 to n-1 using list comprehension, then returns a list of the matrix exp_nbeing raised to the elements of the ascendingly increasing list. i.e.

exp_n**[0, 1, ..., n-1] = [exp_n**0, exp_n**1, ..., exp_n**(n-1)]

So I was wondering if there's a more numpythonic way of doing it(in order to make use of Numpy's broadcasting ability) like:

exp_n**np.arange(1,n,1) = np.array(exp_n**0, exp_n**1, ..., exp_n**(n-1))
6
  • you can do exactly what you described in your last code block. exp_n ** np.arange(1, n) just works. Am I missing something?
    – amdex
    Jul 25, 2019 at 15:06
  • If you want to apply it directly to your matrix, you do need to use broadcasting, as follows: exp_n[None, :, :] ** np.arange(1, n)[:, None, None].
    – amdex
    Jul 25, 2019 at 15:09
  • exp_n ** np.arange(1, n) with exp_n having a dim of (2,2) gives me an error message ValueError: operands could not be broadcast together with shapes (2,2) (4,) though
    – Darren Ng
    Jul 25, 2019 at 15:20
  • You are aware that you go from 1 -> n-1 right? Not from 0 Jul 25, 2019 at 15:37
  • yea. Also your exp_n[None, :, :] ** np.arange(1, n)[:, None, None] works, it's exactly what I wanted. Also it turns out that exp_n[None, :] ** np.arange(1, n)[:, None, None] works too and is a little bit faster. If you put it down as official answer I'll accept it.
    – Darren Ng
    Jul 25, 2019 at 15:38

1 Answer 1

2

You're speaking of a Vandermonde matrix. Numpy has numpy.vander


def c_matrix_vander(n):
    exp = np.exp(1j*np.pi/n)
    exp_n = np.array([[exp, 0], [0, exp.conj()]], dtype=complex)
    return np.vander(exp_n.ravel(), n, increasing=True)[:, 1:].swapaxes(0, 1).reshape(n-1, 2, 2)

Performance

In [184]: %timeit c_matrix_vander(10_000)
849 µs ± 14.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [185]: %timeit c_matrix(10_000)
41.5 ms ± 549 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Validation

>>> np.isclose(c_matrix(10_000), c_matrix_vander(10_000)).all()
True
3
  • Amazing. Thanks.
    – Darren Ng
    Jul 25, 2019 at 15:59
  • In any case, there's actually a exp_n**0 as the first matrix. However the 0 elements in exp_n raising to 0 becomes 1 in numpy, so I'd actually get [[1,1], [1,1]] instead of the identity. I resolve this issue by np.concatenate((II, c_matrix_vander(n)) where II is the identity matrix. Do you see a better way of doing this?
    – Darren Ng
    Jul 25, 2019 at 16:15
  • I'm a bit confused as to what you mean, because 0**0 is 1 mathematically speaking Jul 25, 2019 at 16:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.