Rank of a Potential Evaluation Matrix

In [51]:
import numpy as np
import matplotlib.pyplot as pt

Let's make two particle collections: sources and targets

In [52]:
sources = np.random.randn(2, 200)
targets = np.random.randn(2, 200)

pt.plot(sources[0], sources[1], "go")
pt.plot(targets[0], targets[1], "ro")
Out[52]:
[<matplotlib.lines.Line2D at 0x7f7914197e10>]

Now let's assume each of these points has a charge, and evaluate the potential at each of the other points.

In [53]:
all_distvecs = sources.reshape(2, 1, -1) - targets.reshape(2, -1, 1)
dists = np.sqrt(np.sum(all_distvecs**2, axis=0))
interaction_mat = 1/dists

pt.imshow(dists)
Out[53]:
<matplotlib.image.AxesImage at 0x7f790fe6d198>

Finding the Rank: Attempt 1

How do we find the rank? Get the matrix to echelon form, look for zero rows.

Bonus Q: Is this the same as LU?

In [54]:
from m_echelon import m_echelon
M, U = m_echelon(interaction_mat)
pt.imshow(np.log10(1e-15+np.abs(U)), cmap="gray")
pt.colorbar()
Out[54]:
<matplotlib.colorbar.Colorbar at 0x7f790fd858d0>

Finding the Rank: Attempt 2

In [55]:
U, sigma, V = np.linalg.svd(interaction_mat)

pt.semilogy(sigma)
Out[55]:
[<matplotlib.lines.Line2D at 0x7f790ffdf908>]
In [56]:
k = 60
Uk = U[:, :k]
Vk = V.T[:, :k].T

Ak = (Uk * sigma[:k]).dot(Vk)

np.linalg.norm(interaction_mat - Ak, 2)
Out[56]:
7.0952237471847974
In [56]: