This isn't a complete algorithm, but rather an idea that hopefully someone else can build off of.
One potential way to frame the problem is that we seek two polynomials, $P(z)$ and $Q(z)$:
$$
P(z) = \sum_{i=1}^{N} z^{x_{i}} \qquad
Q(z) = \sum_{i=1}^{N} z^{y_{i}}
$$
Such that their product corresponds to our sums:
$$
R(z) = P(z) \cdot Q(z) = \sum_{i=1}^{N} \sum_{j=1}^{N} z^{s_{i,j}}
$$
You could then factor $R(z)$, using a polynomial time factoring algorithm to obtain its irreducible factors over the integers. Then you'll need to combine these factors into two polynomials with $N$ non-zero terms each.
Note that $x_1 = y_1 = 0$, $s_{1,1} = 0$, meaning all factors will have a constant term, and that there will be no overall constant factor.
Also note that $P(1) = Q(1) = N$. This implies that for each of our polynomial factors $p_i(z)$, we have $p_i(1) \mid N$. This can allow us to efficiently mark factors as mutually exclusive (one belonging to each of $P(z)$ and $Q(z)$) if $p_i(1) \: p_j(1) \nmid N$. Polynomials can then be multiplied efficiently using the FFT, and terms checked. In many cases, there won't even be that many polynomial factors- almost all polynomials with coefficients zero or one are irreducible over the integers.
I'm struggling to bound this step in polynomial time, but maybe someone else can weigh in and prove a bound on the number of irreducible factors a polynomial with coefficients of ones and zeros can have.
One important caveat: the time complexity of this approach will be dependent on the values of your sequence, not just their length. For large $s_{i,j}$ this algorithm will perform terribly.