Accelerating scientific python program

Question

I have the following code in python:

def P(z, u0):
    x = np.inner(z, u0)
    tmp = x*u0
    return (z - tmp)


def powerA2(A, u0):
    x0 = np.random.rand(len(A))
    for i in range(ITERATIONS):
        x0 = P(np.dot(A, x0), u0)
        x0 = x0 / np.linalg.norm(x0)
    return (np.inner(np.dot(A, x0), x0))

np is numpy package.

I am interested in running this code for matrices in size of 100,000 * 100,000, but it seems that there is no chance for this program to run fast (I need to run it many times, about 10,000).

Is there any chance that tricks like multi-threading would work here?

Does anything else help to accelerate it?

Answer 1

You could consider using Pythran. Compiling the following code (norm.py):

#pythran export powerA2(float [][], float[])
import numpy as np

def P(z, u0):
    x = np.inner(z, u0)
    tmp = x*u0
    return (z - tmp)

def norm(x):
    return np.sqrt(np.sum(np.abs(x)**2))

def powerA2(A, u0):
    ITERATIONS = 100
    x0 = np.random.random(len(A))
    for i in range(ITERATIONS):
        x0 = P(np.dot(A, x0), u0)
        x0 = x0 / norm(x0)
    return (np.inner(np.dot(A, x0), x0))

with:

pythran norm.py

yields the following speedup:

$ python -m timeit -s 'import numpy as np; A = np.random.rand(100, 100); B = np.random.random(100); import norm' 'norm.powerA2(A, B)'
100 loops, best of 3: 3.1 msec per loop
$ pythran norm.py -O3 -march=native
$ python -m timeit -s 'import numpy as np; A = np.random.rand(100, 100); B = np.random.random(100); import norm' 'norm.powerA2(A, B)'
1000 loops, best of 3: 937 usec per loop