Scalar fields can be propagated through non-paraxial systems using the Gaussian beam decomposition method. However, for high NA objectives, this scalar treatment is not sufficient to correctly describe the electromagnetic fields inside the focal region due to high ray bendings, which result in a significant change in the polarization state of light. To model these vectorial effects, the Gaussian beam decomposition method has to be extended to include the polarization state of light. In this work we have combined it with the three dimensional polarization ray tracing in order to propagate vectorial fields through high NA optical systems. During the Gaussian beam decomposition, the polarization state of each individual beamlet is represented by a polarization vector [𝐸𝑥, 𝐸𝑦, 𝐸𝑧 ] associated with its central ray. Individual Gaussian beams are then propagated through the system using the complex ray tracing method. The effect of the optical system on the polarization state of each beam is computed by applying the three dimensional polarization ray tracing of the corresponding central rays. Finally the individual beams are superposed coherently in the plane of interest resulting in the complete vectorial field. We apply the proposed method to compute the vectorial field inside the focal region of a high NA microscope objective lens and compare our result to the vectorial Debye integral method. We find that the Gaussian beam decomposition method overcomes serious limitations of algorithms relying on Fourier transforms, i.e. the field sampling requirements are less critical in high NA focusing and in the presence of large aberrations. However, sharp edges in the amplitude profile are difficult to handle as they should be replaced with smooth Gaussian edge.