Smooth Pycnophylactic Interpolation is a disaggregation where z_i have minimal quadratic slopes subject to the Pycnophylactic Principle.

Mathematically: minimize $f(z) := (D_x z)^T \cdot (D_x z) + (D_y z)^T \cdot (D_y z)$, subject to: $\forall r: \sum\limits_{i} z_i \cdot q_i^r = Z_r$, where $D_x$ and $D_y$ are the linear operations that result in the partial discrete difference in the $x$ and $y$ directions respectively, i.e.: ${(D_x z)}_{x,y}$

= $z_{x,y} - z_{x-1,y}$ , and $q_i^r$ is the incidence matrix with the amount of raster cell $i$ in region $r$, i.e. $\forall i: \sum\limits_{r} q_i^r = 1$

See also these comments of ChatGPT with a linear algebra solution and an iterative solution. Note that the derivation of $y_{ij}$ on the raster boundary seems incorrect there, but without affecting the final result, and the adaption of the lambdas may be implemented better.

Since

$$(\textbf{D}z) := z_{i} − z_{i−1}$$

one can derive that $$( ((D z)^T D z) = (z^T D^T D z) = \sum\limits_{i1}^n(z_i - z_{i-1})^2\sum\limits_{i1}^n(z_i^2 + z_{i-1}^2 - 2 z_i z_{i-1}) z_0^2 + \sum\limits_{i1}^n(2 z_i^2 - 2 z_i z_{i-1}) - z_n^2 )$$

and $$\frac{\partial f(z)}{\partial z_i} = (z^T D^T D)^T_i + (D^T D z)_i = (D^T D^{TT} z)_i + (D^T D z)_i = 2(D^T D z)_i$$

this convex optimization problem can be reformulated as setting to zero of: $$\frac{\partial[ f(z) + \sum\limits_{r} \lambda_{r} (\sum\limits_{i} z_i *q_i^r - Z_r)]}{\partial z_i}$$

$$= \sum\limits_{i \in { x, y } } 4 z_i - 2 z_{i-1} - 2 z_{i+1} + \sum\limits_{r} \lambda_{r} q_i^r ) ),$$

subject to $$\forall r: \sum\limits_{i} z_i * q_i^r = Z_r$$

from which follows that $$z_{x,y} = {1 \over 4} (z_{x-1,y} + z_{x+1,y} + z_{x,y-1} + z_{x,y+1} ) - {1 \over 8} \sum\limits_{r} \lambda_{r} q_i^r$$

Tobblers work, 1979