Since you are dealing with a homogeneous population and regressing out ancestry components, this may not be as large of an issue as you suspect. It would be nice to see some simulations of the false positive rate you expect when doing this. I imagine there are also a lot of important, true interactions within-chromosome.

Another idea for a way to partition PGSs is by the sign of the PGS SNP effect size. Perhaps if a SNP significantly interacts with negative-effect-size PGS SNPs (but not positive), or vice-versa, this could help with placing causal SNPs/genes relative to other known players in a pathway or distinguishing between a pathways' functions in contrasting diseases. However this may not work since networks are complex and involve maybe activating and inhibitory interactions.

Perhaps you could reduce complexity even further by building a PGS per chromosome, testing all chromosome PGS x chromosome PGS interactions, and then doing a subsequent SNPxSNP GWAS between the SNPs on the implicated chromosomes. I'm not sure how many SNPs are on each chromosome, but this could potentially reduce computational needs and the multiple testing burden since the number of epistatic interactions per trait is very small.

If you didn't regress out ancestry, could the PGS term already sufficiently account for population structure? It may not be necessary to remove ancestry components if the PGS term absorbs polygenic background and therefore ancestry, allowing you to use a larger and more diverse population to estimate the SNPxPGS term. It's unclear to me whether it would fully account for this, but may be something to try.