This paper presents an iterative and decoupled perturbative stochastic Galerkin (SG) method for the variability analysis of stochastic linear circuits with a large number of uncertain parameters. State-of-the-art implementations of polynomial chaos expansion and SG projection produce a large deterministic circuit that is fully coupled, thus becoming cumbersome to implement and inefficient to solve when the number of random parameters is large. In a perturbative approach, component variability is interpreted as a perturbation of its nominal value. The relaxation of the resulting equations and the application of a SG method lead to a decoupled system of equations, corresponding to a modified equivalent circuit in which each stochastic component is replaced by the nominal element equipped with a parallel current source accounting for the effect of variability. The solution of the perturbation problem is carried out in an iterative manner by suitably updating the equivalent current sources by means of Jacobi- or Gauss-Seidel strategies, until convergence is reached. A sparse implementation allows avoiding the refinement of negligible coefficients, yielding further efficiency improvement. Moreover, for time-invariant circuits, the iterations are effectively performed in post-processing after characterizing the circuit in time or frequency domain by means of a limited number of simulations. Several application examples are used to illustrate the proposed technique and highlight its performance and computational advantages.