A preconditioned adaptive integral method is presented for the expedient extraction of the broadband scattering parameters of large arrays of through-silicon vias. The method relies on a proper decomposition of the original kernel functions in the Foldy-Lax based integral equation statement of the pertinent electromagnetic boundary value problem to generate translation-invariant expressions. Thus, Fast-Fourier Transform based convolution can be used for the efficient calculation of the matrix vector products in the Krylov projection-based iterative solution of the sparse matrix equation resulting from the numerical approximation of the governing integral equation. A multi-frontal preconditioner is constructed and implemented for the Krylov projection-based iterative solution of the sparse matrix. The accuracy of the proposed methodology is validated through comparison with HFSS™-obtained scattering parameters of a simple through-silicon via arrays, while its capacity in handling large arrays is demonstrated through its application to the analysis of a 25×25 TSV array.