Checking Rank Estimates¶

In [1]:
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt

Let's make two particle collections: sources and targets

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

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

What's our convergence ratio $\rho$: $$ \rho=\frac{d(c,\text{furthest target})}{d(c,\text{closest source})}? $$

Well, what should $c$ be, really? (Technically, we can use any $c$, so we're free to choose whichever we think is best.)

In [3]:
c = np.sum(targets, axis=1) / targets.shape[1]

c
Out[3]:
array([ 9.98878098, 10.02648225])
In [4]:
pt.plot(sources[0], sources[1], "go")
pt.plot(targets[0], targets[1], "ro")
pt.plot(c[0], c[1], "bo")
Out[4]:
[<matplotlib.lines.Line2D at 0x7f6364f514d0>]

Assemble the interaction matrix¶

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

Estimate the rank depending on precision $\varepsilon$¶

First, obtain an idea of what $\rho$ is:

In [6]:
dist_tgt_to_c = np.sqrt(np.sum((c.reshape(2, 1) - targets)**2, axis=0))
dist_src_to_c = np.sqrt(np.sum((c.reshape(2, 1) - sources)**2, axis=0))

rho = np.max(dist_tgt_to_c) / np.min(dist_src_to_c)
rho
Out[6]:
0.2537137489568833

Then plot the numerical rank depending on epsilon:

In [7]:
_, sigma, V = la.svd(interaction_mat)
pt.semilogy(sigma)
Out[7]:
[<matplotlib.lines.Line2D at 0x7f6364ea8310>]
In [9]:
eps_values = 10**(-np.linspace(1, 12))

For the given precisions $\varepsilon$, find the associated numerical rank:

In [10]:
def numrank(eps):
    return np.sum(sigma > eps)

ranks = [numrank(e) for e in eps_values]
pt.semilogx(eps_values, ranks)
Out[10]:
[<matplotlib.lines.Line2D at 0x7f63659fced0>]

Now compare with our estimate:

In [11]:
pt.semilogx(eps_values, ranks)
pt.semilogx(eps_values, (np.log(eps_values)/np.log(rho)-1)**2)
Out[11]:
[<matplotlib.lines.Line2D at 0x7f6364bb14d0>]
  • We estimated that the rank would grow quadratically.
  • Comments on how good our estimate is?
In [ ]: