Yes, Helmholtz decomposition is inherently a non-local process because it is only unambiguous when the integrals are performed over all of a simply connected space that contains all relevant sources of the field (where a "source" is defined as any region in the field where the divergence or curl is non-zero). As the usage of the magnetic scalar potential demonstrates, it fails on sub-spaces because of the ambiguity inherent in the classification. In standard Helmholtz decomposition, divergenceless and solenoidal are taken as synonymous, likewise with divergent and irrotational (curl free). Fields that have both zero divergence and zero curl are assumed to be excluded by the boundary condition at infinity that the fields vanish there.
When you limit your considerations to finite sized regions, though, the assumption breaks down because a source that is outside of your bounded region can produce a field that is both divergenceless and irrotational everywhere in the region of interest. For that reason, local categorization is only unambiguous if it is done via positive properties. In other words, a vector field can be locally described unambiguously by breaking it down into three components:
\begin{align}
\mathbf{F}(\mathbf{x}) & = \mathbf{F}_{\mathrm{div}}(\mathbf{x}) + \mathbf{F}_{\mathrm{curl}}(\mathbf{x}) + \mathbf{F}_{\mathrm{har}}(\mathbf{x}),\ \mathrm{where} \\
\nabla\cdot\mathbf{F}_{\mathrm{div}}(\mathbf{x}) &\neq 0\ \quad\mathrm{somewhere}, \\
\nabla\times\mathbf{F}_{\mathrm{curl}}(\mathbf{x}) &\neq 0\ \quad\mathrm{somewhere}, \\
\nabla\cdot \mathbf{F}_{\mathrm{har}}(\mathbf{x}) & = 0\ \quad\mathrm{everywhere},\ \mathrm{and} \\
\nabla\times \mathbf{F}_{\mathrm{har}}(\mathbf{x}) & = 0\ \quad\mathrm{everywhere}.
\end{align}
The subscript of $\mathbf{F}_{\mathrm{har}}(\mathbf{x})$ is meant to be short for "harmonic" in analogy to harmonic functions because it can always be expressed in the region of interest as the negative gradient of a harmonic function. You could also call the harmonic term $\mathbf{F}_{\mathrm{ext}}(\mathbf{x})$, since it can also be modeled as being produced by sources external to the region of interest (when that region has finite size).
In this categorization, the reconstruction would work like this:
\begin{align}
\mathbf{F}_{\mathrm{div}}(\mathbf{x}) & = \frac{1}{4\pi} \int_V \frac{\mathbf{x}-\mathbf{x}'}{\left|\mathbf{x}-\mathbf{x}'\right|^3} \nabla'\cdot \mathbf{F}(\mathbf{x}') \operatorname{d}V' \\
\mathbf{F}_{\mathrm{curl}}(\mathbf{x}) & = \frac{1}{4\pi} \int_V \nabla\times\frac{\nabla'\times \mathbf{F}(\mathbf{x}')}{|\mathbf{x}-\mathbf{x}'|} \operatorname{d}V' \\
\mathbf{F}_{\mathrm{har}}(\mathbf{x}) & = \mathbf{F}(\mathbf{x}) - \mathbf{F}_{\mathrm{div}}(\mathbf{x}) - \mathbf{F}_{\mathrm{curl}}(\mathbf{x}).
\end{align}
What this exercise makes clear is if you try to take the limit as $V\rightarrow 0$, the region will either have some kind of delta function source that leaves the field ill defined there (e.g. point charge, line charge, or sheet charge) or the integrals that define the divergent and solenoidal components vanish, leaving only a locally harmonic field. Thus, the classification is only useful when done in regions of finite size, and only when the region of interest is the entire possible space.
That is why the Wikipedia formulae quoted in the question do not have an ambiguity in their derivation, they have an unstated assumption that no sources beyond the boundary contribute to the local field. As an example, consider a uniform field in the region of interest. Was it produced by one (or more) infinite charged sheet(s) (divergent source), or infinite current sheet(s) (solenoidal source)? Assuming the region of interest is a cube with faces aligned parallel and perpendicular to the field, the Wikipedia sources would have a combination of charged sheets on the faces the field pierces, and a current sheet circulating around the other four.
Point being, the process is only unambiguous when sources are only allowed in a closed region of interest.