Basic linear algebra states that $\det(I-S)$ must be non-zero so that a solution to your linear system exists. On the other hand, if your determinant is (numerically) zero, the basis vectors in your coefficient matrix $(I-S)$ do not span the whole vector space, so that you can't construct a general right-hand side $b$. Think of this picture:

In the example, the rhs $b$ is a 3D-vector with a $z$-component, but the vectors in the coefficient matrix only span the $x-y$-plane. So, the best you can do is to reproduce the vector in the $x-y$-plane, and forget about the $z$-direction. That is, you essentially solve for the projection $b^{||}$ of the vector $b$ onto the $x-y$-plane. By this you determine only two of the three parameters $x_1$ and $x_2$ in your solution vector. The third parameter $x_3$ can't be determined and is usually set to zero (such that $||x||$ is as small as possible).
To accomplish this, there are several regularization methods available. The two basic ones are the SVD and ridge regression:
The Singular Value Decomposisiton gives an easy method to form the pseudoinverse. It is used to decompose your matrix as $(I-S) = U D V^T$, where $D$ is the diagonal matrix of singular values. Now the pseudoinverse is formed by inverting all entries in $D$, and setting those $D_i$ larger than $1/\epsilon$ to zero (where $\epsilon$ is a small cutoff parameter).
In ridge regression, aka L2- or Tikhonov-regularization, you add a small term $\alpha$ to the diagonal of the moment matrix $(I-S)^T(I-S)$, which then becomes regular so that the linear system can be solved.
In the limit $\alpha,\epsilon \rightarrow 0$, both methods give the same result, namely the solution in the subspace that is spanned by the vectors in $(I-S)$.
By the way, to clarify: this is no solution, but an approximation to the solution. So the answer of @FredericiPoloni still holds. This here is just a workaround.