Ambiguity in applying Newton's shell theorem in an infinite homogeneous universe Newton's shell theorem has two corollaries:


*

*The gravitational attraction of a spherically symmetric body acts as if all its mass were concentrated at the center, and

*The gravitational acceleration inside the cavity of a spherically symmetric body vanishes.
Consider a spaceship floating freely in space. In a homogeneous universe, the combined attraction from all matter should cancel out, and the spaceship should stay motionless. Nevertheless, I'm free to divide the attraction into several parts originating from different parts of the universe: In the figure below, I've divided the universe into a red sphere centered on some arbitrary point (×) with my spaceship located at the edge of the sphere, plus infinitely many shells centered on the same point.
By corollary #1, the gravitational attraction of the red sphere equals that of all its mass centered at the point ×. By corollary #2, the combined acceleration of the spaceship from all mass in the green shell vanishes. The same can be said for the blue shell, the orange shell, and so on ad infinitum.
Hence my spaceship should start accelerating toward ×. By choosing the sphere large enough, I should be able to make it accelerate arbitrarily fast, and by choosing the location of × I can make it accelerate in any direction.

Of course this doesn't work, but why?
My best guess is that, even in an infinite universe, you can't keep adding spheres because you'll exit the observable universe, in which case there's no way to feel the gravity in part of the shell so that it's no longer symmetric. Perhaps also the expansion of the universe matters. But see the last two points below.
A few more things to consider:


*

*The mass of the red sphere increases with the chosen radius $r$ as $r^3$, while the acceleration it generates is proportional to $r^{-2}$; hence the acceleration increases linearly with the chosen $r$.

*Our universe — the "Universe" — has an average density of some $10^{-29}\,\mathrm{g}\,\mathrm{cm}^{-3}$. Hence if I set $r$ equal to the radius of the observable Universe (46.3 billion light-years), the acceleration is a minuscule $10^{-7}\,\mathrm{cm}\,\mathrm{s}^{-2}$. If that bothers you, choose another universe where $\rho$ is ten orders of magnitude higher.

*Our Universe is not really homogeneous, but on large enough scales ($\gtrsim$ half a billion light-years) it seems it is. Still, the acceleration of the spaceship will be dominated by nearby sources. If that bothers you, choose a sufficiently homogeneous universe.

*On the scales we're considering, the Universe is not governed by Newtonian dynamics, but by general relativity. If that bothers you, use Birkhoff's theorem instead — I think the issue is the same.

*If the issue is really that the size of the observable universe matters, then my intuition tells me that I can just choose an arbitrarily old universe where the asymmetric contribution from the most distant shells is arbitrarily small.

*If the issue is that the universe expands (so that gravity from the far side of a shell is somehow weakened, or "redshifted"), then my intuition tells me that I can just choose a sufficiently static universe.
 A: Updated 07.11
We can chose the model to discuss the problem and so let us chose:
Model: Newtonian mechanics/Newtonian gravity, with the Universe filled with uniformly dense matter, interacting only gravitationally (in cosmology this called “dust matter”), and at the initial time of our spaceship journey all this matter is at rest.

Hence my spaceship should start accelerating toward ×. By choosing the sphere large enough, I should be able to make it accelerate arbitrarily fast, and by choosing the location of × I can make it accelerate in any direction.

Absolutely!

Of course this doesn't work, but why?.

It does work. If we assume that initially the spaceship was at rest together with the whole universe it will reach the point × in time needed for the ship to fall into a point mass equal to the mass of the pink sphere. 
The problem is that by that time all of the pink sphere also falls toward that same point as well, as do all other colored spheres  and the rest of the universe also. If our astronaut checks its distance to the point × before the spaceship falls into it she would notice that this distance has decreased, but at the same time is she checks her surroundings she would notice that the spaceship is surrounded by precisely the same matter particles that when the journey started only they are closer to each other and to the spaceship. This distance contraction is simply a Newtonian version of Big Crunch event. 
If the universe is filled with matter interacting only gravitationally and we assume that the density of matter will stay uniform throughout the universe, then the only conclusion would be that such universe is not static. It has either (Newtonian version of) Big Bang in its past or Big Crunch in its future (or in our model, since we chose initial moment as a turning point from expansion to contraction, it has both).  
It may seem that the whole Universe falling toward our chosen point × is an absurdity, since we have chosen this point arbitrarily. But in this situation there is no paradox, the acceleration of all matter toward this point is due to the fact that in our setup there is no “absolute space”, no set of outside stationary inertial observers which could give us  absolute accelerations, instead we can only choose a reference point × (or rather specify an observer located at this point and at rest with respect to surrounding matter) and calculate relative accelerations toward this point. 
Recall, that the first principle of Newtonian mechanics states that every particle continues in its state of rest or uniform motion in a straight
line unless it is acted upon by some exterior force. For an isolated system, for example collection of gravitating objects of finite total mass we could (at least in principle) place an observer at rest so far away that it could be considered an inertial object. This would allow us to define a reference frame with respect to which we would measure accelerations. But in our Newtonian cosmology matter is filling the whole Universe, there is no observer on which gravity is not acting, so there is no set of reference frames defined by observers “at infinity” only observers inside the matter concentrations that are affected by the gravitational forces. 
While there is no absolute accelerations, the relative positions ($\mathbf{d}_{AB}(t)= \mathbf{x}_A(t)-\mathbf{x}_B(t)$ between objects $A$ and $B$ comoving with the matter of the universe) do have a meaning independent of the choice of reference point. This relative positions, relative velocities ($\dot{\mathbf{d}}_{AB}$), relative accelerations, etc. constitute the set of unambiguously defined quantities measurable within our universe.

then my intuition tells me that I can just choose a sufficiently static universe.

This intuition is wrong, if there is a gravitational force that would accelerate your spaceship toward ×, then it would also be acting on a nearby matter (call them dust particles or planets or stars) producing the same acceleration, so all of the universe would be falling toward ×.

Note on Newtonian cosmology it may seems that Newtonian theory of gravitation is ill suited to handle homogeneous spatially infinite distributions of matter. But one can try to separate the physics of the situation from the deficiencies of particular formalism and possibly to overcome them. As a motivation we could note that over large, cosmological distances our universe to a high degree of accuracy could be considered spatially flat, and the velocities of most massive objects relative to each other and to the frame of CMB are very small compared with the speed of light, meaning that Newtonian approximation may be appropriate. While we do know that general relativity provides a better description for the gravitation, Newtonian gravity is computationally and conceptually much simpler. This seems to suggest that it is worthwhile to “fix” whatever problems one encounters while attempting to formalize cosmological solutions of Newtonian gravity. 
The most natural approach is to “geometrize” Newtonian gravity and instead of “force” consider it a part of geometry, dynamical connection representing gravity and inertia. This is done within the framework of Newton–Cartan theory.
As a more detailed reference, with an emphasis on cosmology, see this paper (knowledge of general relativity is required):


*

*Rüede, C., & Straumann, N. (1996). On Newton-Cartan cosmology. arXiv:gr-qc/9604054.


Newton–Cartan theory underscores conceptual similarities between Newtonian gravity and general relativity, with Galilei group replacing the Lorentz group of GR. 
The general approach is coordinate-free and is closely related to the machinery of general relativity, but a specific choice of local Galilei coordinates would produce the usual equations for acceleration ($\mathop{\mathrm{div}} \mathbf{g} = - 4\pi \rho$), with gravitational acceleration now being part of Newtonian connection. Homogeneous and isotropic cosmological solutions are a straightforward lifts of FLRW cosmologies.
While equations are the same, we may already answer some conceptual questions. 


*

*Since gravitational acceleration is part of the connection, there is no reason to expect it to be an “absolute” object, there would be gauge transformations that would alter it. We can have multiple charts on which we define the physics with the normally defined transition maps between. 

*We can have a closed FRW cosmology, the “space” does not has to be a Euclidean space, it could be torus $T_3$ (field equations require that locally the space is flat). Since the spatial volume of a closed universe varies, and tend to zero as the universe approaches the Big Crunch, this asserts that not just matter but space itself collapses during the Big Crunch (to answer one of the comments). 

*It is quite simple to include the cosmological constant / dark energy thus making the models more realistic.
Note on answer by user105620: If we formulate a regularization procedure by introducing a window function $W(\epsilon,x_0)$ that would make potential well behaved. This provides us with an another way to “fix” problems of our cosmological model. The acceleration of our spaceship computed with this regularization is indeed dependent on the choice of $x_0$ in the limit $\epsilon\to 0$, which is the consequence of the same freedom in choosing the reference point ×. But he/she just should not have stopped there. Divergences   requiring the use of regulators and  ambiguities remaining after regularization are quite normal features in developing physical models. The next step would be identifying the physically meaningful quantities and checking that those are independent on the regulator artifacts. In our case neither potential $\Phi$ nor gravitational acceleration $\mathbf{g}$ are directly observable in this model. Relative positions, relative velocities and relative accelerations are  observable and those are turning to be independent of the regulator parameter $x_0$. 
A: 
by choosing the location of × I can make it accelerate in any direction.

This freedom of choice is the key to the puzzle. I'll assume Newtonian gravity in a static universe filled with a homogeneous dust. 
Let the ship be at the origin. The ship feels a force proportional to $x$ towards the centre of the sphere of radius $x$ centered at $\pmb{x}$, but it also feels the exact opposite force towards the centre of the identical but disjoint sphere centered at $\pmb{-x}$, so these two forces cancel exactly. In each case, I'm only considering the mass inside the ball and ignoring the mass outside it, per the shell theorem.
The same logic applies to any arbitrary $\pmb{x}$.
A: I'd like to address, in a rigorous way, what's going on mathematically that leads to this apparent contradiction. Newton's shell theorem, as proved by Newton, is a statement about the gravitational field as defined through the Newton's law of universal gravitation,
$$\mathbf{g}(\vec{x}) = \int_{\mathbb{R}^3}\rho(\vec x') \frac{(\vec x'-\vec x)}{|\vec x'-\vec x|^3} d^3x'. \tag{1}$$
Where $\rho: \mathbb{R}^3 \to \mathbb{R}_+$ is the mass-density function, which we will take to be constant. Whether this formula is formally what one wants to call Newtonian Gravity or not, this is where our contradiction must lie. By definition, the above formula implies that the $i$th component $\mathbf{g}_i(\vec x)$ of the gravitational field is 
$$\mathbf{g}_i(\vec x) = \rho \int_{\mathbb{R}^3} \frac{x_i'-x_i}{|\vec x'-\vec x|^3} d^3x',$$
and now our integrand is simply a real-valued function, a situation with which we are comfortable. However, the fundamental issue with this expression is that, though it looks like we may well call it zero by symmetry, the integrand is not integrable in the Lebesgue or improper Riemann sense because it is not absolutely integrable, i.e.
$$\int_{\mathbb{R}^3} \frac{|x_i'-x_i|}{|\vec x'-\vec x|^3} d^3x' = \infty$$
in the Legesgue sense. Here's the kicker: because our integrand is not integrable, we cannot expect theorems indicating consistency under change of coordinates and passing to iterated integrals to apply. But this is precisely our issue: each time you apply the shell theorem about a different choice of center, you are invoking a change to a particular set of spherical coordinates and computing the resulting expression via an iterated integral (one must, as Newton's shell theorem applies to an "infinitesimally" thin spherical shell). Because of the above technical issue, the values obtained in each case need not be consistent with each other.
As discussed by user105620, different types of issues arise in the formulation of Newtonian gravity through a potential, wherein $\mathbf{g}$ is determined by the conditions $\vec \nabla \cdot \mathbf{g} = \rho$, $\vec \nabla \times \mathbf{g} = 0$, and a boundary condition on $\mathbf{g}$. If $\rho$ does not decay sufficiently fast (as in the hypotheses of the linked result), this formulation is not generally well-posed, i.e. such a $\mathbf{g}$ may not exist (though, if it does, it is probably unique, depending on the boundary condition). 
Existence aside, the shell theorem in this case, proven by the divergence theorem, hinges on being able to assume spherical symmetry of $\mathbf{g}$ from that of $\rho$. One can easily show that this works fine for the standard case of $\rho$ decaying sufficiently fast with the boundary condition $\mathbf{g} \to 0$ at infinity, but it is not clear at all how to prescribe a physically reasonable boundary condition that ensures it is allowed otherwise. Indeed, for the constant $\rho$ case, $\mathbf{g}(\vec x) = \frac{\rho}{3} (\vec x - \vec x_0)$ satisfies the PDE conditions for any $\vec x_0$, but such solutions do not differ by a constant, so the linked uniqueness statement above implies that all standard types of boundary conditions (Dirichlet, Neumann, and mixed) can only select one of these. That is, in potential Newtonian gravity, standard choices of boundary conditions cannot generically allow us to assume spherical symmetry of $\mathbf{g}$ from that of $\rho$ when $\rho$ does not decay, and hence the shell theorem generally fails in this case.
Ultimately, then, your contradiction comes down this: considering the two most basic theories of Newtonian gravity which naturally include the shell theorem, it turns out that one theory simply doesn't make mathematical sense in the non-decaying $\rho$ case, while the other theory's shell theorem necessarily breaks down in the non-decaying $\rho$ case.
A: From a very quick skim it seems the existing answers are excellent, so I'll instead contribute some of the physics and philosophy literature. I too was concerned by this issue after reading a certain paper (Peacock 2001, incidentally), until I discovered centuries of thought preceded me!
Your concern was apparently first raised by Bishop Berkeley, in discussion with Newton himself. Much later Seeliger (1890s) sharpened and popularised the critique. See Norton (1999), "The cosmological woes of Newtonian gravitation theory" for history. Norton also discusses the analogous issue for Coulomb's law of electric force.
Remarkably, Newtonian cosmology was only worked out after the general relativistic case, by Milne and also McCrea. Here I particularly mean the rate of expansion, which closely resembles the relativistic Friedmann equations incidentally. [I'm assuming a homogeneous and isotropic universe. Otherwise, see Buchert & Ehlers (1997).] But again your objection was raised. Finally, Heckmann & Schucking (1955) are credited with making Newtonian cosmology great again rigorous.
Norton was yet another who independently raised the centuries-old objections. Malament (1995) defended by describing 3 formulations of Newtonian gravity: the $1/r^2$ force law, Poisson's equation, and Newton-Cartan theory. Norton (1995) concurred, yet added that acceleration becomes relative! Tipler (1996a, 1996b) has nice papers from the same time. Wallace (2017) looks interesting, such as the section title "2. Non-uniqueness of solutions to Poisson's equation".
A: The problem lies in the boundary conditions. Ignoring factors of $G$ and $\pi$, gauss's law of gravitation relates the gravitational potential $\Phi$ to the mass density $\rho$ by
$$\rho=-\nabla^2 \Phi. $$
In order to have a unique, well-defined solution, we need to specify boundary conditions for $\Phi$. Usually, we assume that $\rho$ dies off sufficiently quickly at spatial infinity that a reasonable choice of boundary condition is $\Phi(|\vec x|\to\infty)=0$ is. The shell theorem relies on this assumption. However in your example $\rho$ does not die off at infinity and is instead non-zero everywhere and therefore the shell theorem fails.
Often when a given scenario in physics doesn't, but almost, satisfies the 'if' part of a theorem, it can be helpful to  try and modify the problem so that it does. Therefore we can use a window function $W_\epsilon(x-x_0)$ that dies off quickly as $x\to\infty$ but $\lim_{\epsilon\to0} W_\epsilon =1$ to regulate the charge density. [e.g. take $W_\epsilon(x-x_0)=e^{-\epsilon (\vec x-\vec x_0)^2}$.] Then we can replace your uniform charge density $\rho$ by
$$\rho\to\rho_{\epsilon,x_0}\equiv \rho W_\epsilon(x-x_0) .$$
In this case, the shell theorem does hold. However, the result we get is not regulator-independent, that is if we solve for $\Phi_{\epsilon,x_0}$ using the charge distribution $\rho_{\epsilon,x_0}$ and then send $\epsilon \to0$, we find that our answer still depends on choice of $x_0$. 
This is the mathematically rigorous way to see that there really is an ambiguity when applying the shell theorem to such a situation! 
Edit: There seems to be some debate in the comments as to whether the shell theorem should be proved with forces or with Gauss's law. In reality, it doesn't matter, but I will address what goes wrong if you just use forces. Essentially, Newton's laws are only guaranteed to be valid if there is a finite amount of matter in the universe. Clearly if there is uniform mass density throughout all of space, then there is an infinite amount of matter, so the shell theorem fails. The requirement that $\rho(|\vec x|\to \infty)\to 0$ 'sufficiently quickly' from above is more precisely that $\int d^3 x \rho(x) <\infty$, which is just the condition that there is a finite amount of matter in the universe. 
A: You should sum up contribution of shells centered at the test particle, not at the center of the arbitrary sphere. The shell theorem is not applicable in this case.
Let's consider a simpler 1D case.

Suppose a rope pulling (tug of war) game with infinite number of players in the both teams with a test particle at zero.
All players produce the same pull and are located at all integer distances from zero, but at zero itself and at $1$ the players are missing.
Which team will win? Here we step into the area of divergent series and integrals. The cardinalities of the sets {-1,-2,-3,...} and {2,3,4,...} are the same, but they have different numerocities. By definition, numerocity is the sum of the indicator function over the set: $\sum_{k=-\infty}^\infty p(k)$, where $p(k)$ is $1$ if $k$ is in the set and $0$ otherwise.
So, the numerocity of the left team is $\sum_{k=-\infty}^{-1} 1$ while the numerocity of the right team is $\sum_{k=2}^\infty 1$. The both sums are divergent, but their regularized values are different.
The regularized value of the left-hand sum is $-1/2$ while the regularized value of the right-hand sum is $-3/2$.
They also can be expressed with integrals as $\int_{-\infty}^{-1/2}dx$ and $\int_{3/2}^\infty dx$.
If the infinite parts cancel each other, the finite parts give victory to the left-hand-side team.
In a more complicated 3D case one should account for situations, where the infinite parts may not cancel each other (this is not the case of your original question, but may arise if at the sides of a test particle are  infinite cones, cylinders, paraboloids, sponges or other infinite shapes, which need to be compared).
Due to Fourier or Laplace transform we can formally denote the divergent integral $\int_0^\infty dx$ as $\pi\delta(0)$ or, for simplicity, $\tau$.
This way, the strength of the left-hand team in our example is $\int_{-\infty}^{-1/2}=\tau-1/2$ and the strength of the right-hand team is $\int_{2/2}^\infty=\tau-3/2$
The integrals of monomials, needed in multi-dimensional complicated cases will follow the following law:
$\int_0^\infty x^n dx=\frac{\left(\tau +\frac{1}{2}\right)^{n+2}-\left(\tau -\frac{1}{2}\right)^{n+2}}{(n+1)(n+2)}$
For more information on divergent integrals, see here (need MathML-capable browser (Firefox, PaleMoon).
A: It is the case that by invoking this red circle with center x, to use the shell theorem you must be physically at the edge of the universe for it to work. Because if this is not the case, you arbitrarily chose a circle and point x on the left side of your spaceship, which you could have equally as well have drawn to the right side of you in this homogenous universe. And those cancel out, similar to the trivial proof of Newton in corollary 2 (proposition 70). For his proof uses actual spheres, with actual centers; not arbitrary ones. That is the problem with your scenario.
A: The spaceship doesn't accelerate in either direction but the gravitational effect between points can be thought of as slowing down the expansion of the universe. There is a lecture here by Susskind using precisely this approach to derive the Friedmann Equations:
https://www.youtube.com/watch?v=938_TNP4aUs
