Decreasing the time necessary to enter the coefficients of a matrix

Question

Using Python, I want to build a square matrix whose coefficients are functions evaluated at some given points. The matrix is lower triangular, yet it takes more than 200 seconds to enter all the coefficients when the size is 1000. This is quite surprising to me as I have already dealt with square matrices of size greater than 1 million.

I guess that the loss of time comes from the function that defines all the non-zero coefficients of the matrix: it is a product of four terms, each of them involving a complex power (see the minimal example below). I would be really happy if someone could help me: the purpose is to reach a size of 1 million within a few seconds (hope this is doable...). Maybe the use of ** is not optimal for the complex power; or maybe there is a quicker way to fill a lower triangular matrix? Any help will be appreciated!

For the code below, Ns is the size of the square matrix.

import math
from math import *
import numpy as np
from numpy import exp, log, sqrt, cos, sin, arange
import time


tmps1 = time.time()


Ns = 100
h = 1/float(Ns)
Ns = int(Ns)

mu = 0.01

rn = -6.34
rc = 0.86
r1 = 1.32
r2 = 4.16

P = np.zeros((Ns + 1, 1))
for j in range(0, Ns + 1):
    P[j] = r1 + j*(r2-r1)*h

kappan = 0.24
kappac = -0.24
kappa1 = 0.095
kappa2 = -0.095


z = complex(0.01,0.01)


def E(exponent,p):
    return (  ( ((p-rn)/(2-rn))**((exponent)/(2*kappan)) )*( ((p-rc)/(2-rc))**((exponent)/(2*kappac)) )*( ((p-r1)/(2-r1))**((exponent)/(2*kappa1)) )*( ((r2-p)/(r2-2))**((exponent)/(2*kappa2)) )  )

def D(p,r):
    return (      ( 1/(2*complex(0.0,1.0)*(z-mu/r1)) )*(   ( 1/(2*(kappa1-complex(0.0,1.0)*(z-mu/r1))) )*( E(2*(kappa1-complex(0.0,1.0)*(z-mu/r1)),r)*E(2*complex(0.0,1.0)*(z-mu/r1),p) )   -   ( 1/(2*kappa1) )*E(2*kappa1,r)   )      )


A = np.zeros((Ns-1, Ns-1), dtype=np.complex_)

for j in range(1, Ns-1):
    for k in range(1, j+1):
        A[j,k] = D(P[j+1], P[k+1]) - D(P[j+1], P[k])


tmps2 = time.time()-tmps1
print "\n\nExecution time = %f\n\n" %tmps2

Answer 1

Optimization, Round-1:

Avoid recomputation of D.

d = np.zeros((Ns+1, Ns+1), dtype=np.complex_)
for j in range(2, Ns):
    for k in range(1, j+1):
        d[j,k] = D(P[j], P[k])

A = np.zeros((Ns-1, Ns-1), dtype=np.complex_)
for j in range(1, Ns-1):
    for k in range(1, j+1):
        A[j,k] = d[j+1, k+1] - d[j+1, k]

Reduces time to half

Optimization, Round-2:

Vectorized the computation of D.
Replace computation of d with:

d = np.zeros((Ns, Ns), dtype=np.complex_)
for shift in range(0, Ns-1):
    x = D(P, np.roll(P,shift))
    for j in range(shift+1, Ns):
        d[j,j-shift] = x[j]

Optimization, Round-3:

Vectorized the computation of A.
Replace the computation of A with:

d = np.roll(d, -1, axis=(0,1))
A = d - np.roll(d, 1, axis=1)
A *= np.tri(*A.shape)

Takes around 1.2 seconds for Ns=1000

Edit:

Round-4:

P = np.fromfunction(lambda j, i: r1 + j*(r2-r1)*h, (Ns+1, 1), dtype=float)
Q = np.fromfunction(lambda j, i: r1 + ((j-i)%(Ns+1))*(r2-r1)*h, (Ns+1, Ns), dtype=float)
# P = Q[:,0].reshape(11,1)
d = D(P, Q)
d = np.fromfunction(lambda r, c: d[r,(r-c)], (Ns,Ns), dtype=int)

d = np.roll(d, -1, axis=(0,1))
A = d - np.roll(d, 1, axis=1)
A = np.tril(A)

Thanks Fabio Lipreri for suggestion of fromfunction, I have exploited it to the hilt!

Time Taken: 400ms, for Ns=1000

Answer 2

In addition to the previous answer of ckedar, you can speed up a little by improving the calculation of P vector using numpy vectorization. In the following code I used numpy's fromfunction function:

const = (r2-r1)*h
P = np.fromfunction(lambda j, i: r1 + j*const, (Ns + 1, 1), dtype=float)

Using this code on my machine, with an Ns = 1000, takes around 15.2 µs to generate P against 421 µs using your code (96% speed up).