The quantity that refuses to change: discovering conservation laws from raw data
The problem
Push a child on a swing and watch. The height changes constantly. The speed changes constantly, and reverses twice a swing. But one particular mix of the two never changes at all — high and slow at the top, low and fast at the bottom, always trading one for the other in exactly the same proportion.
Knowing that mix is worth more than any number of photographs of where the swing was. A photograph tells you where it is. The mix tells you everywhere it can ever be, and everywhere it cannot. That is the difference between recording a system and understanding it, and it is the same difference whether the system is a swing, a chemical reactor, a market, or a patient's vital signs. Something that holds steady while everything around it moves is called an invariant, and finding one is the classical shape of discovering a law.
The usual way to look for one is to guess. Write down a family of formulas you think the answer might belong to, fit their coefficients to your data, and see which fits best. That works when you already know roughly what you are looking for, and fails silently when you do not, because the best fit within the wrong family still looks like an answer.
This walkthrough guesses nothing. You hand it logged measurements and it hands back the combination that never moves, or tells you there isn't one. The second half is the important half.
You need numpy, and for the last step a running HeatherDB.
The idea, on the swing
Suppose you had four hundred instruments pointed at the swing, all badly calibrated and none of them measuring anything sensible. Instrument 87 reads some arbitrary function of position and speed. Instrument 212 reads a different arbitrary one. Nobody designed them and none of them is height or energy or anything with a name.
Now log the swing for a while, and for every moment write down not the readings but how much each reading changed since the moment before. You get a huge table: one row per moment, four hundred columns of changes.
Here is the whole trick. If some mix of the instruments is conserved, then that mix's changes must cancel out at every single moment — that is what conserved means. So you are looking for a way of weighting the four hundred columns such that every row of the table adds up to zero. Find that weighting and you have found the invariant, without ever having said what you thought it would look like.
There is a standard tool for "find the weightings this table cannot move". It takes a table and reports the directions it stretches most, down to the directions it barely touches at all, each with a size. The ones with a size of essentially zero are the mixes the table leaves alone, and those are the candidates.
Two precautions keep this from being a trick on yourself. Subtract each instrument's average first, so that an instrument stuck at a constant value cannot win by being trivially constant. And divide by how much each instrument varies, so the answer is not just "whichever instrument was quietest". Both are one line, and without them the pipeline finds a law every time.
Step 1 — make the data
Three systems, and none of them ever shows the pipeline its rule of motion. A frictionless spring, a pendulum, and the same spring with friction added.
def record(accel, amplitude, n=400, dt=0.05):
"""Let go from rest at `amplitude` and log 400 moments of position and speed."""
q, p, out = amplitude, 0.0, []
for _ in range(n):
p += 0.5*dt*accel(q, p); q += dt*p; p += 0.5*dt*accel(q, p)
out.append((q, p))
return np.array(out)
Seven releases from seven heights are used to look for the law. Four different heights are held back and used only for grading, so a law that merely memorised the recordings it was shown has nowhere to hide. The real energy of each system is computed too, but only ever to mark the answer — it is never an input and never available to the step that chooses between candidates.
Step 2 — look for the mixes nothing moves
changes = np.vstack([measure(r[1:]) - measure(r[:-1]) for r in recordings])
_, sizes, mixes = np.linalg.svd(changes / spread, full_matrices=False)
the sizes the search reports
harmonic spring biggest 34.25 smallest three 9.3e-15 9.1e-15 8.8e-15
pendulum biggest 25.98 smallest three 8.9e-15 8.8e-15 8.5e-15
spring + friction biggest 30.85 smallest three 1.1e-14 1.1e-14 1.1e-14
Read that carefully, because it is the trap. Every system reports mixes with a
size of about 1e-14, which is zero as far as a computer is concerned —
including the one with friction, which has nothing conserved at all. There
are always more ways to weight four hundred instruments than there are things
for the data to push against, so unmoved mixes exist whether or not any of them
means anything.
So the sizes are not the answer. They are a shortlist. The answer has to be graded somewhere else.
Step 3 — grade it on data it never saw
Two numbers per system. How closely the discovered quantity tracks the real energy on the four held-back releases — releases at energies absent from the search entirely. And how much it wanders during a single recording compared to how much it differs between recordings, which should be near zero for something that holds still while still telling one recording from another.
graded on releases the search never saw
harmonic spring matches the real energy at 0.978, wanders 0.000 -> a conserved quantity
pendulum matches the real energy at 1.000, wanders 0.000 -> a conserved quantity
spring + friction matches the real energy at 0.465, wanders 0.000 -> nothing conserved
The pendulum is the one to look at twice. Its conserved quantity is not a simple sum of squares — it involves a cosine — and no straight-line fit of position and speed could ever express it. It comes back at 1.000 from raw logged numbers, at heights that were never in the search, because the four hundred deliberately meaningless instruments were nonlinear, so a straight-line mix of them is a curved function of the swing.
Notice also that "wanders" reads 0.000 for the friction case too. A quantity that is nearly constant everywhere trivially holds still along every recording, including the wrong ones. One statistic can always be satisfied by a degenerate answer, which is why the verdict needs both numbers and why the grading happens on recordings the search never touched.
Step 4 — the control is the result
The third row is not a footnote. It is the reason the first two rows mean anything.
Friction drains the spring, so there is no conserved energy for the pipeline to find, and it must come back empty. It does: 0.465, against 0.978 and 1.000. The pipeline is identical in every respect — same instruments, same search, same candidate selection, same grading — and only the answer differs, in the direction the physics demands.
A procedure that returns a law for every input has not found a law. It has found a way to always produce output, which is worth nothing, and it is the default outcome of anything with enough free parameters. The 1.000 on the pendulum is evidence only because the 0.465 was available and came up. If you build one thing from this post, build the control.
Step 5 — keep the law, and repair a bad measurement of it
A law you found once is worth storing, because measuring it properly is expensive. Write it into a collection, and later hand the memory a hasty measurement to be cleaned up.
call("POST", "/db/pendulum_rig/collections/laws/write",
{"vectors": [[float(x) for x in mix]], "metadata": [{"system": "pendulum"}]})
Two short, noisy recordings are a bad basis for anything — twenty-five to forty moments each with measurement noise on top, instead of seven clean recordings of four hundred. Reading the resulting mess against the memory returns whatever stored law it lands nearest.
from seven long clean recordings: 1.000
two 40-step recordings, noise 0.02: 0.152 -> 1.000 after one read
two 30-step recordings, noise 0.05: 0.224 -> 1.000 after one read
two 25-step recordings, noise 0.08: 0.262 -> 1.000 after one read
pure noise, not a law at all: 0.358 -> 1.000 after one read
The repair is real: a measurement worth 0.152 comes back worth 1.000. So is the catch, on the last line. Hand the memory four hundred random numbers that are not a law at all and it returns the pendulum law just as happily, because with one law stored there is nothing else for it to return. The memory repairs; it does not verify. Keep the two jobs apart — cleaning up is a read, and deciding whether the thing you cleaned up is right is the correlation grade from Step 3.
The whole thing
import base64, json, urllib.request
import numpy as np
NF = 400 # how many readings to take of a state
rng = np.random.default_rng(0)
W = rng.standard_normal((NF, 2)) * 2.5 # each reading is a random direction …
B = rng.uniform(0, 2*np.pi, NF) # … with a random offset
SCALE = np.array([1/np.pi, 1/2.0])
def measure(S):
"""Turn a position-and-speed pair into 400 arbitrary readings of it. None of
them knows what system it is looking at; they are deliberately meaningless."""
return np.sqrt(2.0/NF) * np.cos(np.dot(np.atleast_2d(S) * SCALE, W.T) + B)
def record(accel, amplitude, n=400, dt=0.05):
"""Let go from rest at `amplitude` and log 400 moments. The rule of motion is
used here to make the data and is never shown to anything downstream."""
q, p, out = amplitude, 0.0, []
for _ in range(n):
p += 0.5*dt*accel(q, p); q += dt*p; p += 0.5*dt*accel(q, p)
out.append((q, p))
return np.array(out)
def search(recordings):
"""Every step of every recording gives one row: how much each of the 400
readings changed. A quantity that never changes is a mix of readings whose
changes always cancel — so we want the mixes that this table cannot move."""
changes = np.vstack([measure(r[1:]) - measure(r[:-1]) for r in recordings])
levels = np.vstack([measure(r) for r in recordings])
levels -= levels.mean(0) # so a reading that is simply flat can't win
changes -= changes.mean(0)
spread = levels.std(0) + 1e-9 # so a quiet reading can't win by being quiet
_, sizes, mixes = np.linalg.svd(changes / spread, full_matrices=False)
return mixes, sizes, spread
value = lambda mix, S, spread: np.dot(measure(S) / spread, mix)
def drift(mix, recordings, spread):
"""How much the quantity wanders during one recording, against how much it
differs between recordings. Near zero means it holds still but still tells
one recording from another."""
within = np.mean([np.std(value(mix, r, spread)) for r in recordings])
return within / (np.std(value(mix, np.vstack(recordings), spread)) + 1e-9)
SYSTEMS = {
"harmonic spring": (lambda q, p: -q, lambda q, p: 0.5*(p**2 + q**2)),
"pendulum": (lambda q, p: -np.sin(q), lambda q, p: 0.5*p**2 + (1 - np.cos(q))),
"spring + friction":(lambda q, p: -q - 0.15*p, lambda q, p: 0.5*(p**2 + q**2)),
}
FIT = [0.4, 0.7, 1.0, 1.3, 1.9, 2.2, 2.4] # releases used to find the quantity
HELD = [0.55, 1.15, 1.6, 2.05] # releases used only to grade it
def characterise(name):
accel, energy = SYSTEMS[name]
fit = [record(accel, a) for a in FIT]
held = [record(accel, a) for a in HELD]
mixes, sizes, spread = search(fit)
# Among the least-moved mixes, keep the one that holds stillest on recordings
# held back from the search. True energy plays no part in this choice.
mix = min((mixes[-k] for k in range(1, 8)), key=lambda m: drift(m, fit[::2], spread))
S = np.vstack(held)
truth = np.array([energy(q, p) for q, p in S])
return mix, spread, sizes, abs(np.corrcoef(value(mix, S, spread), truth)[0, 1]), \
drift(mix, held, spread)
print("the sizes the search reports")
for name in SYSTEMS:
_, _, sizes, _, _ = characterise(name)
print(" %-18s biggest %6.2f smallest three %s" %
(name, sizes[0], " ".join("%.1e" % s for s in sizes[-3:])))
print("\ngraded on releases the search never saw")
for name in SYSTEMS:
_, _, _, corr, d = characterise(name)
verdict = "a conserved quantity" if corr > 0.9 and d < 0.15 else \
"nothing conserved" if corr < 0.7 else "unclear"
print(" %-18s matches the real energy at %.3f, wanders %.3f -> %s" %
(name, corr, d, verdict))
# ── keep the law, and repair a bad measurement of it ────────────────────────
BASE = "http://localhost:6380"
AUTH = base64.b64encode(b"admin:quickstart-password").decode()
def call(method, path, body=None):
req = urllib.request.Request(BASE + path, method=method,
data=None if body is None else json.dumps(body).encode(),
headers={"Content-Type": "application/json", "Authorization": "Basic " + AUTH})
try:
return json.load(urllib.request.urlopen(req))
except urllib.error.HTTPError as e:
if e.code == 409: # already there from a previous run
return {}
raise
call("POST", "/db", {"name": "pendulum_rig", "dimension": NF})
call("POST", "/db/pendulum_rig/collections", {"name": "laws"})
accel, energy = SYSTEMS["pendulum"]
mix, spread, _, corr, _ = characterise("pendulum")
held = [record(accel, a) for a in HELD]
S = np.vstack(held); truth = np.array([energy(q, p) for q, p in S])
grade = lambda m: abs(np.corrcoef(value(m, S, spread), truth)[0, 1])
print("\nfrom seven long clean recordings: %.3f" % grade(mix))
call("POST", "/db/pendulum_rig/collections/laws/write",
{"vectors": [[float(x) for x in mix]], "metadata": [{"system": "pendulum"}]})
noise = np.random.default_rng(3)
for n, eps in ((40, 0.02), (30, 0.05), (25, 0.08)):
short = [record(accel, a, n=n) + eps * noise.standard_normal((n, 2)) for a in (0.8, 1.6)]
m2, _, s2 = search(short)
bad = min((m2[-k] for k in range(1, 8)), key=lambda m: drift(m, short, s2))
fixed = np.array(call("POST", "/db/pendulum_rig/collections/laws/read",
{"query": [float(x) for x in bad]})["result"])
print(" two %d-step recordings, noise %.2f: %.3f -> %.3f after one read"
% (n, eps, grade(bad), grade(fixed)))
junk = np.random.default_rng(1).standard_normal(NF); junk /= np.linalg.norm(junk)
fixed = np.array(call("POST", "/db/pendulum_rig/collections/laws/read",
{"query": [float(x) for x in junk]})["result"])
print(" pure noise, not a law at all: %.3f -> %.3f after one read"
% (grade(junk), grade(fixed)))
What to change for your own logs
Swap record for your data. The pipeline wants a list of recordings, each
an array of moments, each moment a few numbers. Vitals sampled hourly, a
reactor's temperature and pressure, a portfolio's holdings. Nothing below
record knows or cares what those numbers mean.
More than two numbers per moment. Widen the random projection from two columns to however many you have. Everything else is unchanged. A system with several independent conserved quantities reports several unmoved mixes, and the grading step sorts them out.
Your held-back recordings must be genuinely different. Not a random split of the same recordings — different operating conditions entirely. Different release heights here; different patients, different market regimes, different load levels for you. A split within one condition grades memorisation and passes it.
Build a control you know has no law. The easiest one is your own data with the moments shuffled, which destroys the dynamics and keeps everything else. If the pipeline still reports a law on that, the law is coming from your setup and not from your system.
When to trust the answer
Only with the control run. A number from this pipeline on its own is not evidence. It becomes evidence the moment you can show the same pipeline declining on an input where declining is correct.
Only on held-back conditions. Grading on the recordings the search used tells you the search worked, not that anything was discovered.
Only above the gap. 0.978 and 1.000 against 0.465 is a gap you can see. A result at 0.75 is not a weak law, it is an unclear one, and the honest handling is to widen the data or report nothing.
And "no conserved quantity" is not "no law". The spring with friction obeys a perfectly good rule; it just is not a rule about something staying the same. Absence of an invariant is a specific, checkable claim about a system, not a failure to find anything.
More depth: invariant, null model, and Abstention.