Vectorizing three nested loops with Numpy

Question

I have a complex matrix C with dimensions (r, r) as well as a complex vector of size r. I need to compute a new matrix from C and v following this equation:

where K is also a square matrix of dimensions (r, r). Here is the code to compute K with three loops:

import numpy as np
import matplotlib.pyplot as plt

r = 9

#   Create random matrix
C = np.random.rand(r,r) + np.random.rand(r,r) * 1j
v = np.random.rand(r) + np.random.rand(r) * 1j

#   Original loops
K = np.zeros((r, r))

for m in range(r):
    for n in range(r):
        for i in range(r):
            K[m,n] += np.imag( C[i,m] * np.conj(C[i,n]) * np.sign(np.imag(v[i])) )

plt.figure()
plt.imshow(K)
plt.show()

Removing the loop with i is relatively easy:

#   First optimization
K = np.zeros((r, r))

for m in range(r):
    for n in range(r):
        K[m,n] = np.imag(np.sum(C[:,m] * np.conj(C[:,n]) * np.sign(np.imag(v)) ))

but I am not sure how to proceed to vectorize the two remaining loops. Is it actually possible in this case?

TheCoolDrop · Accepted Answer · 2020-10-29 20:21:10Z

3

I had a lot of these of problems and here is how I usually proceeded to find solutions to writing out vectorized code.

Here is what I have noticed about your summation. Cool conclusion is that you probably do not need vectorization at all, as you can express your whole calculation as a single product of 2D matrics. Here comes...

Lets first define following matrix (sorry for lack of Latex notation, Stackoverflow does not support Mathjax) :

A_{i,j} = c_{i,j}.

B_{i,j} = c_{i,j} * sgn(Im(v_i))

Then you can write your summation as:

k_{m,n} = Im( \sum_{i=1}^{r} c_{i,m} * sgn(Im(v_i)) * c_{i,n}^* ) = Im ( \sum_{i=1}^{r} B_{i,m} * A_{i,n}^* ) = Im( \sum_{i=1}^{r} B_{m,i}^T * A_{i,n}^* )

The expression above inside of Im(.) is the by definition of matrix multiplication equivalent to following :

k_{m,n} = Im( (B^T * A^*)_{m,n} )

Which means that your matrix k can be expressed as product of transpose of matrix B and product of matrix A. In your code the matrix matrix A is assigned already to variable C. So the vectorization could be done as follows:

C = np.random.rand(r,r) + np.random.rand(r,r) * 1j
v = np.random.rand(r) + np.random.rand(r) * 1j
k = np.imag( (C * np.sign(np.imag(v)).T @ np.conj(C) )

And you have avoided both nasty loops and convoluted expressions

answered Oct 29, 2020 at 20:21

TheCoolDrop

1,07611 silver badges19 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

JD80121 Over a year ago

Thank you for the step by step guide, it is very useful!

TheCoolDrop Over a year ago

I am glad to be of assistance. I know how intimidating Numpy can sometimes be, so I try to rationalize each of the steps. One more important thing is, once you vectorized scribble down your thoughts in code comments, since vectorized code is very hard to analyze once you forget what you thought while writting it.

Quang Hoang · Accepted Answer · 2020-10-29 20:17:40Z

2

This looks like matrix multiplication:

out = np.imag((C*np.sign(np.imag(v))[:,None]).T @ np.conj(C))

Or you can use np.einsum:

out = np.imag(np.einsum('im,in,i', C, np.conj(C), np.sign(np.imag(v))))

Verification with your approach:

np.all(np.abs(out-K) < 1e-6)
# True

edited Oct 29, 2020 at 20:17

answered Oct 29, 2020 at 20:13

Quang Hoang

151k11 gold badges64 silver badges86 bronze badges

1 Comment

JD80121 Over a year ago

Awesome, this is exactly what I was looking for!

JD80121 · Accepted Answer · 2020-10-29 20:13:17Z

0

I found something that can work for now. However, one loop remains and since the resulting matrix is symetric, there is still some optimization to be made.

Instead of removing the i loop, I removed the two other ones:

K = np.zeros((r, r), dtype=np.complex128)

for i in range(r):
    K += adjointMatrix(C) @ (np.sign(np.imag(v)) * C)

K = np.imag(K)

with:

def adjointMatrix(X):
    return np.conjugate( np.transpose(X) )

answered Oct 29, 2020 at 20:13

JD80121

6191 gold badge7 silver badges15 bronze badges

Collectives™ on Stack Overflow

Vectorizing three nested loops with Numpy

3 Answers 3

2 Comments

1 Comment

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

2 Comments

1 Comment

Comments

Your Answer

Sign up or log in

Post as a guest

Related