Skip to content

Adding the Mexican hat neighborhood function - #2

Merged
andremsouza merged 2 commits into
andremsouza:masterfrom
Arne49:mexican_hat
Jul 28, 2026
Merged

Adding the Mexican hat neighborhood function#2
andremsouza merged 2 commits into
andremsouza:masterfrom
Arne49:mexican_hat

Conversation

@Arne49

@Arne49 Arne49 commented Jun 1, 2026

Copy link
Copy Markdown
Contributor

Additionaly to the Gaussian and Bubble neighborhood function, now also the Mexican hat neighborhood function implemented.

Arne49 and others added 2 commits June 1, 2026 14:25
@Arne49
Adding the Mexican hat neighborhood function
1f2bf6d
@andremsouza
Make the Mexican hat isotropic and fix the cyclic fold-back
d953993 
The Mexican hat was computed as the outer product of two 1-D Ricker
wavelets, mirroring how _gaussian is written. That reasoning holds only
for the gaussian: exp(-(dx^2+dy^2)/2s^2) factors into exp(-dx^2/2s^2) *
exp(-dy^2/2s^2), so its outer product is exactly the isotropic 2-D
gaussian. The Ricker wavelet is not multiplicatively separable.

The product is therefore not a function of the distance between two
nodes, which is what a neighborhood function must be: cf. sqdist(c, i)
in Kohonen (2013) Eq. (5), and the single "lateral distance" abscissa of
Vrieze (1995) Fig. 3. Concretely, in the diagonal quadrants both 1-D
factors are negative and their product flips positive, putting a
spurious excitatory lobe exactly where the function must inhibit. On a
21x21 grid with sigma=3 the old form gave h = +0.165 at (c+2s, c+2s)
where the isotropic form gives -0.055, and its zero-crossing locus was a
cross rather than the circle r = sqrt(2)*sigma.

Rewrite it as (1 - u) * exp(-u) over u = sqdist(c, i) / (2*sigma^2),
normalized so h_cc = 1. It now crosses zero at r = sqrt(2)*sigma and
bottoms out at -exp(-2) at r = 2*sigma.

Also, since both functions now share one distance computation:

- Fix the cyclic fold-back, which only folded offsets above +length/2
  and left those below -length/2 alone. A winner at row 0 of a 10x10
  toroidal map wrapped correctly (h[9] = 0.6065) while a winner at row 9
  did not (h[0] = 0.0). Replaced with the minimum-image convention,
  which folds both tails. This bug predates this branch and affected
  _gaussian; the new function would have inherited it.
- Reject a non-finite or non-positive sigma instead of dividing by zero.
  The bubble still admits sigma = 0, which selects the winner alone.
- Raise ValueError naming the valid options for an unknown
  neighborhood_function, matching weight_initialization, rather than
  letting a bare KeyError escape.
- Reject mode='batch' with the mexican hat. Kohonen Eq. (8) is a
  weighted mean whose denominator sum_j n_j*h_ji is not sign-definite
  for a signed neighborhood function; on a 12x12 grid 49 of 144
  denominators come out negative, which inverts or explodes the update.

Add tests/test_neighborhood.py, asserting closed-form properties rather
than golden values: isotropy, the zero crossing, the minimum, absence of
a positive lobe beyond the zero crossing, and cyclic wrap symmetry in
both directions.
@andremsouza

andremsouza commented Jul 28, 2026

Copy link
Copy Markdown
Owner

Hi @Arne49, thank you for this, and sorry for how long it took me to get to it.

I've reviewed the implementation and pushed a few commits to your branch instead of sending you
through several rounds of review. The feature is going in. What changed is how the function gets
computed in two dimensions, and I want to explain why, because the reasoning is interesting and the
mistake is an easy one to make.

The separability issue

The implementation builds the 2-D neighborhood as the outer product of two 1-D Ricker wavelets:

ax = (1 - dx**2 / sigma**2) * np.exp(-dx**2 / d)
ay = (1 - dy**2 / sigma**2) * np.exp(-dy**2 / d)
return np.outer(ax, ay)

This follows _gaussian, which was the sensible thing to do, since it's the pattern already in the
file. But that pattern works for the gaussian for a special reason:

exp(-(dx² + dy²) / 2σ²)  ==  exp(-dx² / 2σ²) · exp(-dy² / 2σ²)

The exponential turns a sum into a product, so the outer product of two 1-D gaussians really is the
isotropic 2-D gaussian. That's a property of exp alone. The Ricker wavelet doesn't have it: the
polynomial factor (1 - r²/σ²) refuses to separate.

Why it matters in practice

In the diagonal quadrants both (1 - dx²/σ²) and (1 - dy²/σ²) are negative, so their product
flips positive. That's a sign error across a whole region of the map. Measured on a 21×21 grid with
σ=3, centred at (10, 10):

position outer product isotropic form
diagonal, at (c+2σ, c+2σ) +0.165 −0.055
diagonal, at (c+3σ, c+3σ) +0.008 −0.001
global minimum −0.443 −0.135 (= −e⁻², at r = 2σ)
zero-crossing locus a cross (dx=±σ ∪ dy=±σ) a circle, r = √2·σ

So nodes far out along the diagonals get pulled toward the input with up to 16.5% of the winner's
strength, in the region where a Mexican hat should be inhibiting. Along those diagonals the function
never goes negative at all.

What the sources say

Both references the library is built on require the neighborhood to be a function of one scalar
distance:

  • Kohonen (2013), Eq. (5) writes h_ci(t) = α(t)·exp[−sqdist(c,i) / 2σ²(t)], where sqdist(c,i) is
    "the square of the geometric distance between the nodes c and i in the grid".
  • Vrieze (1995), Fig. 3, which is where the Mexican hat comes from as a lateral-interaction
    function, plots interaction against an abscissa labelled simply "Lateral distance". The text
    describes it as: "When a neuron is firing, near by situated neurons are stimulated, diminishing
    with increasing distance to the firing neuron. From a certain distance on inhibition will take
    place that gradually vanishes when the distance becomes very large."     The coefficient is written
    h_{i i_c} = 1/‖i_c − i‖, and Vrieze notes the grid is assumed to be a metric space.

Every formulation collapses (dx, dy) into a single distance first, then applies the profile. The
outer product never does that, so its value isn't a function of any metric on the grid.

The change

_mexicanhat is now the isotropic 2-D form over u = sqdist(c, i) / (2σ²):

u = self._sqdist(c, sigma) / (2 * sigma * sigma)
return (1.0 - u) * np.exp(-u)

normalized so h(c, c) = 1, consistent with _gaussian. It crosses zero at r = √2·σ and bottoms
out at −e⁻² ≈ −0.135 at r = 2σ.

One consequence worth pointing out: this shifts your 1-D profile slightly too. The 2-D
Laplacian-of-Gaussian uses r²/2σ² where the 1-D Ricker wavelet uses x²/σ², which moves the zero
crossing from σ to √2σ even along the axes. Both conventions are standard, and the 2-D one is what a
2-D grid needs.

Two things that were already broken before your PR

Your function copied the distance computation from _gaussian, which is what I would have done, and
inherited two existing bugs along the way. Both functions now share one code path, so I fixed them
there:

  1. The cyclic fold-back only worked in one direction. dx[dx > shape/2] -= shape folds the positive
    tail and leaves offsets below −shape/2 alone. On a 10×10 toroidal map a winner at row 0 wrapped
    correctly (h[9] = 0.6065) while a winner at row 9 did not (h[0] = 0.0). It now uses the
    minimum-image convention, (d + L/2) % L - L/2, which folds both tails. This had been quietly
    breaking toroidal maps for a while, and your PR is how I found it.
  2. A radius of zero divided by zero. Both functions now reject a non-finite or non-positive radius
    with a clear message. _bubble still accepts σ = 0, where "the winner alone" is a sensible
    reading.

One new restriction

mode='batch' now raises a ValueError when combined with mexicanhat. Kohonen's Eq. (8) batch
update is a weighted mean whose denominator is Σⱼ nⱼ·hⱼᵢ, which is only well defined when h is
non-negative. A signed neighborhood function can drive it to zero or past it. On a 12×12 grid I
counted 49 of 144 denominators coming out negative, which inverts the correction or blows it up.
This is a property of the Mexican hat itself, so it would apply to any correct implementation.
Stepwise modes are unaffected.

Also in the pushed commits

  • tests/test_neighborhood.py, asserting closed-form properties instead of golden values: isotropy,
    the zero crossing at √2σ, the minimum at 2σ, absence of a positive lobe past the zero
    crossing, and cyclic wrap symmetry in both directions. 46 tests, all passing.
  • An unknown neighborhood_function now raises ValueError listing the valid names, rather than
    letting a bare KeyError escape.
  • A ## Neighborhood functions section in the README covering all three, with the Vrieze citation.
  • One incidental ruff format change, a trailing comma in weight_initialization, so the file is
    format-clean.

Could you look over the commits before I merge? And thanks again. This turned up a real bug in the
existing gaussian that nobody had noticed, which is more than the PR set out to do.

@andremsouza
andremsouza merged commit 1829ea5 into andremsouza:master Jul 28, 2026
@andremsouza

andremsouza commented Jul 28, 2026

Copy link
Copy Markdown
Owner

One correction to what I said above: I went ahead and merged rather than waiting, since I didn't
want this sitting open any longer on my account after you'd already waited since June.

Your review is still very welcome, and it isn't too late to act on. If something in the isotropic
version looks wrong to you, say so here and I'll open a follow-up. The Mexican hat is one function
behind a name in a dict with tests around it, so changing it after the fact is cheap.

Next up is a larger PR modernizing the packaging, adding a proper test suite and CI, and fixing a
few other things the review turned up along the way. That's where a second pair of eyes would help
most, if you're interested. I'll link it here when it's open.

Thanks again for the contribution and for the nudge.

This was referenced Jul 28, 2026
Merged
Merged
@Arne49

Arne49 commented Jul 29, 2026 via email

Copy link
Copy Markdown
Contributor Author

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

2 participants

@Arne49 @andremsouza