import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt
import scipy.linalg.interpolative as sli
eps = 1e-7
sources = np.random.rand(2, 200)
targets = np.random.rand(2, 200) + 3
pt.plot(sources[0], sources[1], "go")
pt.plot(targets[0], targets[1], "ro")
pt.xlim([-1, 5])
pt.ylim([-1, 5])
pt.gca().set_aspect("equal")
def interaction_mat(t, s):
all_distvecs = s.reshape(2, 1, -1) - t.reshape(2, -1, 1)
dists = np.sqrt(np.sum(all_distvecs**2, axis=0))
return np.log(dists)
def numerical_rank(A, eps):
_, sigma, _ = la.svd(A)
return np.sum(sigma >= eps)
Check the interaction rank:
numerical_rank(interaction_mat(targets, sources), eps)
Idea:
nproxies = 25
angles = np.linspace(0, 2*np.pi, nproxies)
target_proxies = 3.5 + 1.5 * np.array([np.cos(angles), np.sin(angles)])
pt.plot(sources[0], sources[1], "go")
pt.plot(targets[0], targets[1], "ro")
pt.plot(target_proxies[0], target_proxies[1], "bo")
pt.xlim([-1, 5])
pt.ylim([-1, 5])
pt.gca().set_aspect("equal")
Construct the interaction matrix from the target proxies to the targets as target_proxy_mat
.
A note on terminology: The target_proxies
are near the targets but stand in for far-away sources.
target_proxy_mat = interaction_mat(targets, target_proxies)
Check its numerical rank and shape:
numerical_rank(target_proxy_mat, eps)
target_proxy_mat.shape
Now compute an ID (row or column?):
idx, proj = sli.interp_decomp(target_proxy_mat.T, nproxies)
Find the target skeleton as target_skeleton
, i.e. the indices of the targets from which the remaining values can be recovered:
target_skeleton = idx[:nproxies]
Check that the ID does what is promises:
P = np.hstack([np.eye(nproxies), proj])[:,np.argsort(idx)]
tpm_approx = P.T @ target_proxy_mat[target_skeleton]
la.norm(tpm_approx - target_proxy_mat, 2)
Plot the chosen "skeleton" and the proxies:
pt.plot(sources[0], sources[1], "go")
pt.plot(targets[0], targets[1], "ro", alpha=0.05)
pt.plot(targets[0, target_skeleton], targets[1, target_skeleton], "ro")
pt.plot(target_proxies[0], target_proxies[1], "bo")
pt.xlim([-1, 5])
pt.ylim([-1, 5])
pt.gca().set_aspect("equal")
What does this mean?
Can we come up with an equivalent of a multipole expansion?
Check that this works for 'our' sources:
imat_error = (
P.T.dot(interaction_mat(targets[:, target_skeleton], sources))
-
interaction_mat(targets, sources))
la.norm(imat_error, 2)
nproxies = 25
angles = np.linspace(0, 2*np.pi, nproxies)
source_proxies = 0.5 + 1.5 * np.array([np.cos(angles), np.sin(angles)])
pt.plot(sources[0], sources[1], "go")
pt.plot(targets[0], targets[1], "ro")
pt.plot(source_proxies[0], source_proxies[1], "bo")
pt.xlim([-1, 5])
pt.ylim([-1, 5])
pt.gca().set_aspect("equal")
Construct the interaction matrix from the sources to the source proxies as source_proxy_mat
:
A note on terminology: The source_proxies
are near the sources but stand in for far-away targets.
source_proxy_mat = interaction_mat(source_proxies, sources)
source_proxy_mat.shape
Now compute an ID (row or column?):
idx, proj = sli.interp_decomp(source_proxy_mat, nproxies)
source_skeleton = idx[:nproxies]
P = np.hstack([np.eye(nproxies), proj])[:,np.argsort(idx)]
tsm_approx = source_proxy_mat[:, source_skeleton].dot(P)
la.norm(tsm_approx - source_proxy_mat, 2)
Plot the chosen skeleton as well as the proxies:
pt.plot(sources[0], sources[1], "go", alpha=0.05)
pt.plot(targets[0], targets[1], "ro")
pt.plot(sources[0, source_skeleton], sources[1, source_skeleton], "go")
pt.plot(source_proxies[0], source_proxies[1], "bo")
pt.xlim([-1, 5])
pt.ylim([-1, 5])
pt.gca().set_aspect("equal")
Check that it works for 'our' targets:
imat_error = (
interaction_mat(targets, sources[:, source_skeleton]) @ P
-
interaction_mat(targets, sources))
la.norm(imat_error, 2)