Skeletonization using Proxies¶

Copyright (C) 2020 Andreas Kloeckner

MIT License Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

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

import scipy.linalg.interpolative as sli

def interp_decomp(A, k):
    idx, proj = sli.interp_decomp(A, k)
    P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]
    return P, idx[:k]

eps = 1e-7
In [2]:
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")
No description has been provided for this image
In [3]:
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)
In [4]:
def numerical_rank(A, eps):
    _, sigma, _ = la.svd(A)
    return np.sum(sigma >= eps)

Check the interaction rank:

In [5]:
numerical_rank(interaction_mat(targets, sources), eps)
Out[5]:
np.int64(9)

Idea:

  • Don't want to build whole matrix to find the few rows/columns that actually matter.
  • Introduces "proxies" that stand in for
    • all sources outside the targets or
    • all targets outside these sources

Target Skeletonization¶

In [6]:
nproxies = 25

angles = np.linspace(0, 2*np.pi, nproxies)
target_proxies = 3.5 + 1.5 * np.array([np.cos(angles), np.sin(angles)])
In [7]:
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")
No description has been provided for this image

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.

In [8]:
target_proxy_mat = interaction_mat(targets, target_proxies)

Check its numerical rank and shape:

In [9]:
numerical_rank(target_proxy_mat, eps)
Out[9]:
np.int64(24)
In [10]:
target_proxy_mat.shape
Out[10]:
(200, 25)

Now compute an ID (row or column?):

In [16]:
P, target_skeleton = interp_decomp(target_proxy_mat.T, nproxies)

Check that the ID does what is promises:

In [21]:
tpm_approx = P.T @ target_proxy_mat[target_skeleton]

la.norm(tpm_approx - target_proxy_mat, 2)
Out[21]:
np.float64(4.318300865373527e-15)

Plot the chosen "skeleton" and the proxies:

In [22]:
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")
No description has been provided for this image

What does this mean?

  • We have now got a moral equivalent to a local expansion: The point values at the target skeleton points.
  • Is it a coincidence that the skeleton points sit at the boundary of the target region?
  • How many target proxies should we choose?
  • Can cheaply recompute potential at any target from those few points.
  • Have thus reduce LA-based evaluation cost to same as expansion-based cost.

Can we come up with an equivalent of a multipole expansion?

Check that this works for 'our' sources:

In [25]:
imat_error = (
    P.T @ interaction_mat(targets[:, target_skeleton], sources)
    -
    interaction_mat(targets, sources))

la.norm(imat_error, 2)
Out[25]:
np.float64(5.426832317627278e-09)

Source Skeletonization¶

In [26]:
nproxies = 25

angles = np.linspace(0, 2*np.pi, nproxies)
source_proxies = 0.5 + 1.5 * np.array([np.cos(angles), np.sin(angles)])
In [27]:
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")
No description has been provided for this image

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.

In [28]:
source_proxy_mat = interaction_mat(source_proxies, sources)
In [29]:
source_proxy_mat.shape
Out[29]:
(25, 200)

Now compute an ID (row or column?):

In [31]:
P, source_skeleton = interp_decomp(source_proxy_mat, nproxies)
In [32]:
tsm_approx = source_proxy_mat[:, source_skeleton] @ P

la.norm(tsm_approx - source_proxy_mat, 2)
Out[32]:
np.float64(3.471111818432676e-15)

Plot the chosen skeleton as well as the proxies:

In [33]:
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")
No description has been provided for this image

Check that it works for 'our' targets:

In [37]:
imat_error = (
    interaction_mat(targets, sources[:, source_skeleton]) @ P
    -
    interaction_mat(targets, sources))

la.norm(imat_error, 2)
Out[37]:
np.float64(5.7679856249083905e-09)
  • Sensibly, this is just the transpose of the target skeletonization process.
    • For a given point cluster, the same skeleton can serve for target and source skeletonization!
  • Computationally, starting from your original charges $x$, you accumulate 'new' charges $Px$ at the skeleton points and then only compute the interaction from the source skeleton to the targets.

Hierarchical Skeletonization¶

In [40]:
gathered_skeletons = np.concatenate(
    [targets[:, target_skeleton],
     targets[:, target_skeleton] + np.array([0, 1]).reshape(-1, 1),
     targets[:, target_skeleton] + np.array([1, 0]).reshape(-1, 1),
     targets[:, target_skeleton] + np.array([1, 1]).reshape(-1, 1),
    ], axis=1)

pt.plot(gathered_skeletons[0], gathered_skeletons[1], "ro")
Out[40]:
[<matplotlib.lines.Line2D at 0x7f23a29f3620>]
No description has been provided for this image
In [38]:
target_proxies = 4 + 2 * np.array([np.cos(angles), np.sin(angles)])
pt.plot(target_proxies[0], target_proxies[1], "bo")

parent_proxy_mat = interaction_mat(gathered_skeletons, target_proxies)

idx, proj = sli.interp_decomp(parent_proxy_mat.T, nproxies)
parent_skeleton = idx[:nproxies]

pt.plot(gathered_skeletons[0, parent_skeleton], gathered_skeletons[1, parent_skeleton], "ro")
Out[38]:
[<matplotlib.lines.Line2D at 0x7f23a4a9f380>]
No description has been provided for this image
In [ ]: