The estimation of large sparse inverse covariance matrices is a ubiquitous statistical problem in many application areas such as mathematical finance, geology, health, and many others. The -regularized Gaussian maximum likelihood (ML) method is a common approach for recovering inverse covariance matrices for datasets with a very limited number of samples. A highly efficient ML-based method is the quadratic approximate inverse covariance (QUIC) method. In this work, we build on the advancements of QUIC algorithm by introducing a highly performant sparse version of QUIC (SQUIC) for large-scale applications. The proposed algorithm focuses on exploiting the potential sparsity in three components of the QUIC algorithm, namely, construction sample covariance matrix, matrix factorization, and matrix inversion operations. For each component, we present two approaches and provide supporting numerical results based on a set of synthetic datasets and a stylized financial autoregressive model. Testing conducted on a single modern multicore machine show that using advanced sparse matrix technology, SQUIC can recover large-scale inverse covariance matrices of datasets with up to million random variables within minutes. In comparison to competing ML-based algorithms, SQUIC is orders of magnitude faster with comparable recovery rates.