Assuming that your kernel is somewhat smooth, use low-rank approximation.
Here's a naive example:
import numpy as np
N=2000
input=np.random.random(N)
x=np.linspace(-1,1,N)
y=np.linspace(-2,2,N)
X,Y=np.meshgrid(x,y,sparse=True)
A = np.exp(1j*2*np.pi*X*Y)
output = np.dot(A, input)
U,S,V = np.linalg.svd(A)
# find truncation rank for given tolerance
k = np.nonzero(S < 1e-12)[0][0]
print('Truncation rank:', k)
Xk = U[:,:k].dot(np.diag(S[:k]))
Yk = V[:k,:]
output2 = Xk.dot(Yk.dot(input))
# check quality of approximation
print(np.linalg.norm(output - output2, np.inf))
%timeit np.dot(A, input)
%timeit Xk.dot(Yk.dot(input))
Output:
Truncation rank: 32
6.3535593536e-12
100 loops, best of 3: 3.9 ms per loop
10000 loops, best of 3: 61.4 µs per loop
What the above does is compute a truncated singular value decomposition of your matrix and use the resulting low-rank approximation of A to speed up the matrix-vector product. So, instead of using $A$, you use an approximation
$$
A \approx X_k Y_k^T,
\qquad X_k,Y_k \in \mathbb R^{N\times k}
$$
with $k \ll N$.
In the example, this leads to a speedup of more than 60x, and this will only get better as N gets larger. That's because the full matrix-vector product has complexity $\mathcal O(N^2)$, whereas the low-rank matrix-vector product has complexity $\mathcal O(Nk)$.
The result is not exact, but you can influence the accuracy by playing with the truncation rank.
Of course, computing the SVD is very slow. Instead, there are faster alternatives for computing low-rank approximations to a matrix, for instance: fast randomized SVD or, what I would recommend, Adaptive Cross Approximation (ACA). The latter is a black-box algorithm for computing a low-rank approximation without ever even computing the entire matrix A and thus is very well suited for very large problems. It's also easy to implement.
The original literature on ACA (in particular, by M. Bebendorf from 2000 on) is somewhat inaccessible. Instead, here's a random paper which describes ACA in a basic form: (link).
If your kernel is relatively simple, it might even be possible to derive a separable approximation to it analytically, which directly leads to a low-rank approximation instead of having to compute it algorithmically. Check the literature on "multipole expansion" and "fast multipole method" for many examples of this.
Finally, the bottleneck in your original program was actually computing A, not the matrix-vector product. ACA will eliminate that, since it doesn't need the full matrix.
UPDATE:
I just realized that scipy already includes a fast method for approximating the SVD in the little-known scipy.linalg.interpolative module. Here's an example:
import scipy.linalg.interpolative as sli
U,S,V = sli.svd(A, 1e-12)
print('Rank:', U.shape[1])
Xk = U.dot(np.diag(S))
Yk = V.conj().T
output3 = Xk.dot(Yk.dot(input))
# check quality of approximation
print(np.linalg.norm(output - output3, np.inf))
And the output:
Rank: 37
2.62815545927e-10
The difference in timing between exact SVD and approximate SVD (again, this will get even better with larger N due to the involved asymptotics):
%timeit np.linalg.svd(A)
%timeit sli.svd(A, 1e-12)
1 loop, best of 3: 11.4 s per loop
1 loop, best of 3: 297 ms per loop
However, this still requires you to compute the entire matrix A first and doesn't get rid of that bottleneck. Nevertheless, if you need to compute A only once and then have many matrix-vector products, it may be a suitable option.
