When thinking about abstraction in the physical world, one of the main questions is: to what extent does the world factor into pieces with relatively-low-dimensional interactions between them? In some classes of systems, we can do a quick-and-dirty check pretty easily, and this post shows an example.
Our example will be an n-body potential problem, between two somewhat-separated clusters of particles. We’ll look at the gravitational potential induced by one cluster of particles (at positions (x_1, ..., x_n) \in \mathbb{R}^{3 n}), at each of the positions of the particles in the other cluster ((y_1, ..., y_m) \in \mathbb{R}^{3 m}). The potential is given by the usual formula for a gravitational or electric charge system:
\phi(y_j, x) = - \sum_i \frac{1}{|y_j - x_i|}
Then the question is: can we factor \phi(\_, x) via a low-dimensional summary (at least for the points y of interest), i.e. find some g, h for which
\phi_y(x) := (\phi(y_1, x), ..., \phi(y_m, x)) = (g(h(x), y_1), ..., g(h(x), y_m))
and the codomain of h is \mathbb{R}^k for some k << 3n?
(In this case, the answer is at least partially known: the fast multipole algorithm would calculate the potentials at each y_j via a low-dimensional summary of x. On the other hand, fast multipole is not (as far as I know) known to use the smallest possible summary for a given level of precision, so we could view the method below as a quick check on how close to optimality fast multipole is.)
What answer should we expect? Well, suppose the two clusters of particles are very far apart. Then we should be able to use a point-mass approximation, so a 3-dimensional summary should suffice (i.e. center of mass position of the x-cluster. For clusters only slightly closer together, a point-mass plus dipole should suffice. At the other limit, if the two clusters are overlapping, then we’d expect lots of local interactions between the particles, and therefore possibly no low-dimensional summary at all. The interesting empirical question is: what does the summary dimension look like for intermediate distances? We’ll test it below by looking at one cluster randomly distributed in the unit cube, and the other cluster randomly distributed in the unit cube translated by (2, 2, 2), each with about 1000 points.
Then the quick-and-dirty test is:
- calculate the Jacobian of \phi_y(x)
- take its singular vector decomposition
- check how many singular values are not close to zero
The idea here is that, if g, h exist and are differentiable, then
\frac{\partial \phi_y}{\partial x} = \frac{\partial g}{\partial h} \frac{\partial h}{\partial x}
which is a low-rank decomposition (of rank k) of \frac{\partial \phi_y}{\partial x}. So, using the SVD, we check for a low-rank decomposition.
Here are the resulting singular values for one run:
We get a staircase pattern! This suggests that there are “plateaus” of precision in this problem - we need only a handful of components to get relative precision to within ~e^{-2} (the rough size of the first drop), then accounting for more components mostly doesn’t add precision until we add enough to “fall off” the plateau and down to the next level. By the time we account for ~120 components, the precision “bottoms out” at numerical error; adding more no longer helps at all.
When would this kind of method fail? There will be various tricky cases, but I expect the main two potential problems in practice are:
- Discreteness - obviously a method based on differentiation won’t work for discrete functions
- Chaos - in chaotic systems, derivatives typically diverge. The more useful analogous check would differentiate one distribution with respect to another, but that’s a far more computationally difficult problem. And unfortunately, chaos is probably responsible for rather a lot of the nice abstractability in the real world - think e.g. the high-level state variables of an ideal gas.
Here’s the code:
import numpy as np
from matplotlib import pyplot
n = 500
m = 450
class Sim:
def __init__(s):
s.n = n
s.pts = np.random.random((s.n, 3)) # x, y, z positions of n pts in the unit cube
def displacements(s, x):
m = x.shape[0]
return np.einsum("ij, k -> ikj", s.pts, np.ones(m)) - np.einsum("ij, k -> kij", x, np.ones(s.n))
def distances(s, x):
displacements = s.displacements(x)
return np.sqrt(np.einsum("ikj, ikj -> ik", displacements, displacements))
def phi(s, x):
distances = s.distances(x)
return np.sum(1.0/distances)
def jac_phi(s, x):
# Returns axis order (x_pts, sim_pts, space_dims)
displacements = s.displacements(x)
distances = s.distances(x) # I am annoyed that cache doesn't like numpy arrays
return - np.transpose(np.transpose(displacements, (2, 0, 1)) / distances ** 3, (2, 1, 0))
sim = Sim()
x = np.array([2.0, 2.0, 2.0]) + np.random.random((m, 3))
jac_phi = sim.jac_phi(x)
# Flatten the x-dimensions
flatter_jac_phi = np.reshape(jac_phi, (m, sim.n*3))
# The main show
u, s, v = np.linalg.svd(flatter_jac_phi)
pyplot.plot(np.log(s[:200]))