10⁴ 1 10⁻⁴ 10⁻⁸ 10⁻¹² two components carry the whole 200 × 200 table everything after them is floating-point zero — flat on to σ₂₀₀ singular values of the 200 × 200 coupling matrix (first 24 of 200) · σ₁ = 5586 · σ₂ = 5279 · σ₃ = 1.1e-10

kernel methods · molecular data · 2026 · report

Predicting excitonic couplings

The learning curve had a knee in it. The knee was the finding.

Code
github.com/tanmayhinge/bimolecular-ml
Status
complete: model built, curve measured, the data's own structure found underneath it
Stack
Python, NumPy, SciPy

The question

Given two molecules, A and B, predict the excitonic coupling between them from their geometries alone. The library holds 200 A-monomers and 200 B-monomers, which pair into 40,000 dimers, and every coupling in the table was paid for with a quantum-chemistry calculation. So the question a model is actually asked here is not whether it can fit 40,000 numbers. It is how few of those calculations you need before the model can supply the rest. That is what a learning curve measures, and it is the deliverable.

The curve I got was not the shape I expected, and chasing the reason turned out to be worth more than the model.

The representation

Every monomer in the library has the same 15 atoms, the same composition (C6H5NO3), and — this is the part worth noticing — the same atomic order. That makes the usual Coulomb-matrix machinery unnecessary. Row k already denotes the same atom in every molecule, so there is nothing to sort and no reason to fall back on eigenvalues. The diagonal, 12 Z2.4, is then identical in all 400 molecules and carries no information at all, so it goes. What is left is the 105 off-diagonal entries of the upper triangle, checked against dscribe.

How a pair of molecules becomes one kernel entry A dimer's two monomers each become a 15 by 15 Coulomb matrix. The diagonal is identical in every molecule and carries no information, so it is discarded, leaving the 105 off-diagonal entries of the upper triangle as the descriptor. A Laplacian kernel scores each monomer separately and the two are multiplied. Because L1 distance is additive across the blocks, that product is the same model as an ordinary kernel on the concatenated descriptors. a dimer A B 15 atoms each Coulomb matrix 15 × 15 descriptor 105 entries 105 entries diagonal dropped Laplacian kernel k(A, A′) k(B, B′) product kernel k(A, A′) · k(B, B′) = one kernel on the concatenation
Each monomer becomes 105 numbers. The two are scored separately and multiplied — which, because L1 distance is additive across the two blocks, is exactly the same model as one kernel on the concatenated descriptors. The natural baseline and the product kernel are not two things.

Model and cost

A pair is scored with the product kernel K((A,B),(A′,B′)) = k(A,A′) · k(B,B′), with k Laplacian, exp(−d/σ), and d the L1 distance between descriptors. On a full grid of training dimers the kernel matrix is a Kronecker product, K = KAKB, so kernel ridge regression can be solved from the eigendecompositions of the two 200 × 200 blocks rather than from the product — O(n3) instead of O(n6).

That is not a micro-optimisation. Fitting all 40,000 dimers takes 8 ms. The dense 40,000 × 40,000 solve it replaces would need a 12.8 GB kernel matrix before it did any arithmetic. Hyperparameters come from cross-validation grouped by monomer, so no monomer leaks across folds, and targets are centred before solving.

Two regimes, and they disagree

There are two different questions hiding in "how much data do you need", and they have to go on the same axis to be compared. Hold out dimers, and you are asking the model to predict a new pairing of molecules it has already met. Hold out monomers, and you are asking it to predict a molecule it has never seen. The horizontal axis is training pairs in both cases, because a pair is what costs a calculation; the 400 monomer geometries are cheap by comparison.

Learning curve for the two regimes, on log-log axes Two curves against the number of training pairs. Predicting unseen pairs of monomers the model has already seen falls from an MAE of 9.4 to 0.008, dropping two orders of magnitude either side of 800 training pairs and passing through the reference band. Predicting pairs of entirely new monomers starts near 22 and only reaches 11 with six hundred times more data, against a target standard deviation of 38. Bands are one standard deviation over five draws. 0.01 0.1 1 10 16 125 500 2,000 8,000 training pairs (reference calculations) mean absolute error reference MAE quoted in the brief standard deviation of the target, 38.4 800 free parameters unseen monomers — new molecules unseen pairs of known monomers
σ and λ are fixed at the selected values rather than re-tuned at each point, so the curve is a statement about training-set size and nothing else. Five random draws per point; bands are one standard deviation. Held-out test error throughout.

The lower curve does what a learning curve should. It passes through the reference band the brief quotes as achievable and ends at an MAE of 0.0083 with 8,000 training pairs. That also settles a question the raw data leaves open: the couplings have a standard deviation of 38.4 and a range of about ±110, so a reference MAE of 0.01 to 0.1 is either quoted in different units or the target is simply very learnable. It is the second. There is no unit mismatch.

The upper curve does nothing at all. Going from 16 training pairs to 10,000 — six hundred times the reference data — moves the error from about 22 to about 11, against a target standard deviation of 38. Same kernel, same descriptor, same solver. It survives a wide search over σ and λ, both kernel families, and standardised descriptors, and it barely moves as the training monomers go from 50 to 180, so it is not a tuning failure. The cause is visible directly in the descriptors: take any monomer's nearest neighbour in Coulomb-matrix space, and the two most similar molecules in the library still differ in their latent coupling factors by about three quarters of the typical magnitude. Two hundred molecules scattered through a 105-dimensional space have no local neighbourhood to interpolate across.

So the 40,000 labels carry far less information than their count suggests. They pin down the 400 monomers in the library and say very little about a 401st. More pairs cannot fix that. Only more molecules could.

The kernel width has no interior optimum

Worth stating plainly, because it makes the fixed σ a deliberate choice rather than a loose end. Widen the kernel and the error keeps falling, while the kernel matrix becomes numerically catastrophic — and the error does not care.

Test MAE against kernel width at 4,000 training pairs. There is no minimum to find.
σ (× median L1)test MAEcond(K)
40.241713.5 × 105
160.060628.3 × 106
640.015351.5 × 108
2560.003862.4 × 109
10240.000973.9 × 1010
40960.000256.3 × 1011

A kernel that keeps improving as it flattens towards a linear model is telling you the target is much simpler than the representation. Re-tuning σ at every point of the learning curve would only track the edge of whatever grid it was given, and would make each point a different model, so the curves above fix σ at 64 × the median L1 distance and λ at 10−10.

The knee

Now the shape. A learning curve in this field is normally a straight line on log-log axes — error falling as a power of training-set size — and the exponent is the number you quote. This one is not straight. It is nearly flat out to a few hundred pairs, drops two orders of magnitude between 500 and 2,000, and flattens again. Fitting a power law would produce a number, and the number would be a summary of something that is not a power law.

The drop lands at about 800. And 800 is 2 × (200 + 200): exactly the number of free parameters if every coupling were a product of two numbers per monomer. That is a hypothesis read off a curve, which is worth very little until it is checked directly — so, checked directly. The coupling matrix is exactly rank 2. Its third singular value is 1.1 × 10−10 against a leading 5,586; that is the cliff in the figure at the top of this page. Fit the best rank-2 matrix agreeing with the observed entries and read off the rest, using no chemistry whatsoever — not the Coulomb matrix, not the kernel, not σ — and every held-out coupling comes back to a maximum error of 10−11.

Which is the real explanation of the knee. The error collapses exactly where enough dimers have been seen to fix 800 numbers.

Putting the two together

The obvious reading — that the chemistry was never needed — is wrong, and the way it is wrong is the useful part. Completion has to identify 800 numbers, and below roughly that many observations the problem is underdetermined: some monomers appear in too few measured pairs to place at all, and the method does not degrade so much as fail. Kernel ridge never fails that way. It has a smooth prior to fall back on, so it is graceful when data is scarce, even though it plateaus early when data is plentiful. The two have precisely opposite strengths, and the crossover sits where the count of observations passes the count of unknowns.

So: seed the completion with the kernel's prediction instead of with zeros. Where observations are plentiful they pin the matrix down and the seed is irrelevant. Where they are scarce and the matrix is underdetermined, the seed decides which of the many consistent completions gets chosen, and the kernel's smooth guess is a good one.

Kernel ridge, rank-2 completion, and the hybrid of the two Three curves against training pairs on log-log axes. Kernel ridge degrades gracefully when data is scarce but plateaus near an MAE of 0.008. Rank-2 completion is useless below roughly 800 observations and exact above it. The hybrid, which seeds the completion with the kernel's prediction, follows whichever parent is better and beats both in the crossover region: at 2,000 pairs it reaches 0.004 where the kernel sits at 0.041. 10⁻¹² 10⁻¹⁰ 10⁻⁸ 10⁻⁶ 10⁻⁴ 10⁻² 1 10¹ 125 500 2,000 8,000 training pairs (reference calculations) mean absolute error 800 free parameters machine precision rank-2 completion (no chemistry) kernel ridge hybrid
The hybrid inherits the better parent everywhere and beats both in the crossover region. Note the vertical axis spans twelve decades; the completion and hybrid curves reach 2 × 10−12, which is floating-point noise, not a small error.
Mean absolute error over three draws. The hybrid never loses to either parent.
training pairskernel ridgerank-2 completionhybrid
1259.3731.29.37
2505.6631.85.65
5001.7033.21.65
1,0000.19510.30.122
2,0000.0410.8570.0036
4,0000.01574.1 × 10−42.4 × 10−7
8,0000.00832.1 × 10−121.9 × 10−12

The crossover region is exactly the one that matters if the point is to spend fewer reference calculations. At 2,000 training pairs the kernel alone sits at 0.041 and the completion at 0.86, while the hybrid is at 0.0036 — an order of magnitude below the better of its two parents. At 4,000 it is at 2.4 × 10−7 while the kernel is still near 0.016. That is the same accuracy from far fewer calculations, which is what putting structure and chemistry together buys over choosing one of them.

The Gaussian process variant

With the same product kernel, a GP and kernel ridge give the same prediction — the posterior mean is the ridge solution, with the noise variance playing the role of λ. What the GP adds is a posterior variance, and the useful question is whether that variance is honest about the split the learning curve just exposed. It should be small for a pair of monomers seen elsewhere in training and large for a pair involving a genuinely new molecule. It is:

Mean posterior variance, trained on a 100 × 100 monomer grid.
test pairs drawn fromposterior variancereading
monomers seen in training1.0 × 10−12effectively certain
monomers never seen3.5 × 10−5seven orders of magnitude less sure

A model that is going to fail on unseen molecules should at least know that it will. This one does, which is the difference between a limitation and a hazard.

Scope

Reproducibility

The notebook runs top to bottom from a clean kernel and is self-checking: the correctness conditions are asserts, not printed output. The Coulomb matrices are checked against dscribe, the Kronecker solve against a dense solve, and the product kernel against the concatenated-descriptor kernel it claims to be identical to. Every figure on this page regenerates from one script against the original data.