This week was the most review-heavy and numerically focused part of the GWPCA work so far.

The branch was no longer at the stage where broad implementation decisions were the main problem. The central question had become:

Can I show, in a way that will stand up to review, that the Python implementation is genuinely aligned with the R reference where it matters?

And just as importantly:

When it is not aligned, can I explain exactly why?

That made this week both technically intense and surprisingly educational.

What Changed This Week

The biggest difference from earlier weeks was that the work stopped being mainly about adding new decomposition functionality and became much more about evidence:

  • evidence that the estimator agrees with GWmodel::gwpca in R on directly comparable outputs
  • evidence that edge cases are handled deliberately rather than accidentally
  • evidence that the package API around GWPCA is stable enough to trust
  • evidence that review comments have been understood and implemented in a way that improves the branch, not just silences a complaint

This is a very different mindset from “feature development,” and I felt that shift strongly while working through the PR.

Concrete Results from the R Comparison

One thing I do not want to lose in a more reflective writeup is the actual observed numbers. So before getting into the interpretation, here are the concrete reference-style comparisons that guided my confidence.

Fixed-bandwidth reference case

For the fixed synthetic reference setup, the local explained-variance ratios matched the R fixture exactly at the precision we were testing.

Explained Variance Ratio (Local PV, %)

Location GWmodel PC1 spml PC1 GWmodel PC2 spml PC2
2 66.79944007 66.79944007 28.21802351 28.21802351
3 63.64566799 63.64566799 31.44977619 31.44977619
4 61.39397786 61.39397786 33.10290676 33.10290676
5 62.47439538 62.47439538 32.70199980 32.70199980
6 57.79939210 57.79939210 32.34985527 32.34985527
7 68.57077309 68.57077309 23.21854269 23.21854269

That already told me something important: the Python implementation was preserving the same local balance of variance across components in the places where the local decomposition was informative.

For the loadings, I compared absolute values because PCA signs are arbitrary. Again, the local structure lined up exactly in the tested reference case.

Component Loadings (Absolute Values)

Location Feature GWmodel PC1 spml PC1 GWmodel PC2 spml PC2
2 A 0.68895623 0.68895623 0.08114459 0.08114459
2 B 0.70852055 0.70852055 0.13404747 0.13404747
2 C 0.15276758 0.15276758 0.98764712 0.98764712
3 A 0.79204176 0.79204176 0.05003965 0.05003965
3 B 0.60487671 0.60487671 0.19889535 0.19889535
3 C 0.08242585 0.08242585 0.97874239 0.97874239
4 A 0.88485219 0.88485219 0.02065589 0.02065589
4 B 0.45039504 0.45039504 0.29329253 0.29329253
4 C 0.11908364 0.11908364 0.95579957 0.95579957
5 A 0.96219121 0.96219121 0.10256587 0.10256587
5 B 0.21097300 0.21097300 0.30528858 0.30528858
5 C 0.17227437 0.17227437 0.94672019 0.94672019
6 A 0.93663008 0.93663008 0.31174872 0.31174872
6 B 0.08511507 0.08511507 0.64500411 0.64500411
6 C 0.33982279 0.33982279 0.69769795 0.69769795
7 A 0.51692673 0.51692673 0.83387714 0.83387714
7 B 0.62819202 0.62819202 0.21598468 0.21598468
7 C 0.58151659 0.58151659 0.50793655 0.50793655

Adaptive-bandwidth reference case

I also did not want the comparison to be accidentally overfit to one neighborhood definition. So I checked an adaptive nearest-neighbor reference case as well.

For the adaptive fixture, the local explained-variance ratios again lined up cleanly:

Location GWmodel PC1 spml PC1 GWmodel PC2 spml PC2
0 67.31107772 67.31107772 24.42745284 24.42745284
1 66.37226587 66.37226587 25.31665071 25.31665071
2 64.99594573 64.99594573 26.88100518 26.88100518
3 63.29265622 63.29265622 30.29365663 30.29365663
4 61.89382622 61.89382622 31.88004401 31.88004401
5 59.07276182 59.07276182 33.25337976 33.25337976

That mattered because fixed and adaptive kernels stress different parts of the implementation. Matching in both cases gave me much more confidence that I was not only lucky in one narrow setup.

R Comparison: What We Actually Observed

The most important thing I want to record from this week is that the R comparison was not just a vague “looks close” exercise. It became much more concrete.

When comparing against GWmodel::gwpca, I observed three different categories of behavior:

1. Things that matched cleanly

These were the encouraging parts:

  • local component directions matched once you account for ordinary PCA sign indeterminacy
  • local explained-variance ratios matched closely
  • the localized structure of the decomposition was telling the same story in both implementations
  • fixed-bandwidth and adaptive-bandwidth comparisons were both encouraging, which mattered because those two modes exercise different neighborhood definitions
  • the locations where the local structure was strongest or weakest were also lining up, so the spatial interpretation was not drifting even when the raw numbers needed closer inspection

This mattered the most scientifically, because those are the outputs users usually interpret first.

Another important practical observation was that the R comparison was not only about one “pretty” location where everything looks good. Once I started checking multiple local neighborhoods, the same pattern kept showing up:

  • the first component was emphasizing the same dominant local structure
  • the second component was capturing the same residual direction
  • the local proportion of variance explained by the first components stayed close enough that the spatial story remained stable

That kind of consistency gave me much more confidence than a single matching table ever could.

2. Things that looked different at first, but for understandable reasons

The biggest example here was the raw eigenvalue scale.

At first glance, a difference in local eigenvalues can feel scary because eigenvalues are central PCA outputs. But once I dug into the comparison more carefully, it became clear that not every mismatch means “the decomposition is wrong.” Some differences are more likely to come from:

  • covariance normalization conventions
  • weighting details
  • whether an implementation reports a quantity before or after a scaling step
  • small numerical conventions around local covariance construction

I also had to keep reminding myself that R-vs-Python comparison in PCA is easy to misread if you are not careful. A few things can create the illusion of disagreement even when the decomposition is effectively telling the same story:

  • one implementation may orient an eigenvector positively while the other orients it negatively
  • local neighborhoods that are nearly degenerate can magnify tiny numerical differences
  • adaptive neighborhoods can differ in subtle ways if tie handling or effective bandwidth scaling is not interpreted the same way
  • raw local eigenvalues are often a less stable cross-language target than local explained-variance ratios

That is why I gradually stopped treating “every number should match exactly” as the only useful success criterion.

That was a very useful lesson for me. Numerical comparison is not just about spotting differences. It is about classifying them correctly.

3. Things that needed to be formalized into tests

This was the real engineering step.

It was not enough to say:

I checked it against R and it seems okay.

Instead, the branch needed tests that would make the comparison durable. So the work shifted into:

  • generating stable reference values
  • deciding which outputs are legitimately comparable
  • handling PCA sign ambiguity cleanly
  • making sure the test strategy itself would survive CI and review

This forced me to define the comparison more carefully than I had at the beginning. In practice, the most defensible direct comparisons became:

  • absolute local loadings, because sign flips are arbitrary
  • local explained-variance ratios, because they are easier to compare meaningfully across implementations than raw eigenvalue magnitude alone
  • both fixed and adaptive neighborhood settings, so the comparison was not accidentally overfit to one mode of operation

I found that once I framed the comparison this way, the branch felt much more scientifically honest.

A More Honest Version of “Parity”

One thing I became more careful about this week was the word parity.

It is tempting to write “perfect parity with R” once you see matching tables for loadings and variance ratios. But the truth is slightly more nuanced, and I think that nuance is important.

What I observed was:

  • strong agreement on the main interpretable local outputs
  • no evidence that the local decomposition structure was fundamentally off
  • some implementation-level details that needed careful framing instead of overclaiming exact equivalence in every raw quantity

That pushed me toward a more defensible standard:

match the R implementation where the outputs are directly comparable, and write tests that make those comparisons explicit.

I think that is a healthier and more scientific way to present the work.

Turning R Comparison into Real Tests

This week also included a lot of practical work around making the comparison testable inside the repository.

The broad shape of that effort was:

  • build fixed and adaptive reference cases
  • compare local loadings and explained-variance ratios
  • treat sign flips correctly instead of misclassifying them as errors
  • keep the fixtures maintainable enough that they do not become a burden

The fixed-vs-adaptive split mattered a lot here.

For the fixed-bandwidth case, I was mainly checking:

  • whether the local weighted covariance construction produced the same effective local structure
  • whether the bisquare weighting behavior was comparable
  • whether the same locations became informative or uninformative under the chosen distance threshold

For the adaptive-bandwidth case, the focus shifted slightly:

  • whether each focal observation was effectively using the intended local neighborhood
  • whether per-location scaling of the neighborhood radius behaved as expected
  • whether local results remained aligned once nearest-neighbor based neighborhoods were involved instead of a fixed spatial radius

This made the comparison feel much more complete than simply saying “I ran one R script and one Python test.”

The branch eventually moved toward embedding the reference values directly in the Python test module rather than making runtime tests depend on external local R artifacts. That ended up being a much better fit for CI and for keeping the PR self-contained.

I liked that outcome because it kept the spirit of the R comparison without making the test suite fragile or dependent on a side-channel setup.

There was also a very practical reason for this change: if the only evidence of R comparison lives in local scripts and external files, then the PR reviewer has to trust that evidence indirectly. Once the relevant reference values are encoded into the Python tests themselves, the evidence becomes part of the branch instead of background context.

How the Tests Were Actually Set Up

This was one of the most important parts of the week, and I do not think the earlier version of this post captured it clearly enough.

I did not want the R comparison to end up as:

  • one local R script
  • two JSON files on my machine
  • and a vague sentence in the PR saying “I checked it”

That would have been hard to review, hard to rerun, and easy to lose.

So the test setup evolved into something much more explicit inside spml/tests/test_decomposition.py.

1. The reference cases were embedded directly into the Python tests

Instead of making CI read external R-generated files at runtime, I moved the relevant reference values into the Python test module itself.

That meant the branch carried the evidence directly:

  • one fixed-bandwidth reference fixture
  • one adaptive-bandwidth reference fixture
  • each including the local loadings and local variance proportions I actually wanted to compare

This was a much better design for the PR because:

  • the tests became self-contained
  • CI no longer depended on local .json or .R files
  • reviewers could inspect the expected values in the same place as the assertions
  • the comparison stopped being “trust me, I ran R earlier”

2. The test data path was made reproducible

I also made the Python side reconstruct the standardized dataset used for comparison in a consistent way, rather than depending on whatever happened to be in memory when I last ran a notebook.

So the setup was roughly:

  • load the reference data
  • standardize it in the same way expected by the GWPCA workflow
  • fit the Python estimator under the same fixed or adaptive neighborhood settings
  • compare only the outputs that are genuinely comparable

That sounds obvious, but it mattered a lot. Cross-language comparison becomes noisy very quickly if the data preprocessing is even slightly inconsistent.

3. I had to normalize the R fixture format before asserting anything

One small but important practical issue was that the R-side reference values were not in a shape that could just be dropped into NumPy assertions.

In particular:

  • some locations were intentionally uninformative and effectively NaN
  • the variance-proportion data needed conversion into a consistent numeric array
  • the loading tables had to be interpreted location-by-location and component-by-component

So part of the Week 6 work was not just “write tests,” but also:

  • write the fixture-conversion helpers
  • make the expected arrays line up with the estimator output layout
  • ensure the comparison logic reflected PCA sign ambiguity instead of misclassifying it as failure

That was a very software-engineering kind of task: not glamorous, but necessary for the numerical story to become reliable.

4. The comparisons were intentionally selective

I did not try to assert exact equality on every raw intermediate quantity. That would have been misleading.

The tests focused on outputs I could defend:

  • absolute local component loadings
  • local explained-variance ratios / local variance proportions
  • both fixed-bandwidth and adaptive-bandwidth cases

And the assertions were chosen for a reason:

  • absolute values for loadings, because PCA sign flips are arbitrary
  • local variance proportions, because they are more stable cross-language targets than raw eigenvalue magnitude alone
  • multiple locations, not just one “nice” focal point

That made the tests much stronger than a single golden-value check.

5. The R-comparison tests were only one layer

Once those reference tests were in place, I also extended the decomposition coverage around the implementation itself.

This ended up including tests for:

  • fit_transform returning the same in-sample values exposed by scores_
  • fit_transform(..., cv=True) correctly forwarding fit-time parameters
  • neighborhoods with zero or degenerate weights
  • mixtures of valid and invalid local fits
  • condition numbers computed from the full local spectrum rather than only the retained components
  • winning_variable_ being missing for failed fits instead of silently wrong
  • unsupervised bandwidth search working without a target y
  • stationarity_test() requiring a fitted estimator and preserving model configuration during permutation refits

That part was important to me because I did not want Week 6 to become only about “R parity.” If a method agrees with R on one fixture but breaks under realistic package usage, that is not enough.

6. CI and review shaped the final setup too

The tests themselves also had to survive the normal project tooling, not just my laptop.

So part of the week was spent making sure the branch held up under:

  • ruff check
  • ruff format
  • type checking
  • the project pytest suite
  • PR review comments about keeping the estimator API and outputs readable

In practice, that meant the numerical tests and the API cleanup were not separate tracks. They kept affecting each other.

For example:

  • once components_ became a pandas-facing object, the tests had to validate the public shape and labeling rather than only raw NumPy contents
  • once decomposition was treated more explicitly as unsupervised, the search and validation tests had to reflect that y is intentionally absent
  • once the R fixtures were inlined, the tests became easier to run in CI and easier to reason about in review

So Week 6 was not only “compare with R.” It was also:

take the R comparison, turn it into maintainable tests, and make those tests fit the package instead of living beside it.

That ended up being one of the most valuable parts of the whole PR.

Reviewer Feedback Was a Big Part of This Week

This was also the week where the review conversation around the PR mattered a lot.

The feedback was not only about numerical comparison. It also touched:

  • use of np.cov instead of unnecessary custom covariance code
  • making decomposition more obviously unsupervised
  • keeping shared validation/search logic generic
  • exposing outputs as readable pandas objects
  • adding tests that check actual values instead of only shapes

Working through those comments was honestly one of the best parts of the week. It forced me to sharpen not just the code, but the reasons behind the code.

There is a big difference between:

  • “I fixed this because the review said so”

and

  • “I now understand why this design is better for the package.”

I felt more of the second by the end of the week.

Edge Cases and Confidence

Once the R comparison work became more solid, I also spent time thinking about confidence in a broader sense:

  • what happens when neighborhoods are degenerate?
  • what happens when some local fits are valid and others are not?
  • what happens when transform-time inputs are misaligned?
  • what happens when stationarity testing refits models under permutations?

That line of thinking eventually fed into stronger edge-case coverage around transform behavior, state reset across repeated fits, and stationarity-test configuration preservation. These are not “R parity” issues directly, but they are part of making the estimator trustworthy rather than only numerically plausible.

This was actually one of the subtle lessons of the week: R parity alone is not enough. A method can match a reference implementation on one carefully prepared dataset and still be brittle in package use. I wanted the branch to do better than that.

Human Side of the Week

This week felt a little like moving from implementation confidence to scientific humility.

I started the project wanting to “match the R implementation.” By this point, I had a more mature understanding of what that really means:

  • know which quantities should match exactly
  • know which ones may differ for understandable reasons
  • know how to test the important parts cleanly
  • know how to communicate the result without exaggeration

That is a much better place to be than simply claiming success too early.

Main Outcome

By the end of Week 6, I felt the branch had crossed an important threshold:

  • the numerical comparison with R was no longer informal
  • the strongest comparisons had been made explicit in tests
  • reviewer concerns had improved the implementation rather than just slowed it down
  • the PR was increasingly about final polish, confidence, and merge-readiness

This was probably the week where the project felt most like real scientific software engineering: not just writing a method, but earning the right to trust it.

<- Back to index