arXiv:2603.02207v3 Announce Type: replace Abstract: The computation of the Log-determinant of large, sparse, symmetric positive definite (SPD) matrices is essential in many scientific computational fields such as numerical linear algebra and machine learning. In low dimensions, Cholesky is preferred, but in high dimensions, its computation may be prohibitive due to memory limitation. To circumvent this, Krylov subspace techniques have proven to be efficient but may be computationally expensive due to the required orthogonalization processes. In this paper, we introduce a Newton-L\'eja interpolation based approach for estimating log-determinants, where the implementation is only based on matrix multiplications and a rough estimation of eigenvalue bounds of the matrix. The method employs an adaptive stopping criterion to select the degree of interpolation according to the data and introduces a centered, scaled-based formulation for computing the logarithmic divided differences. By combining Newton-L\'eja interpolation with a randomized algorithm called Hutch++, we obtain a matrix-free log-determinant estimation framework that avoids explicit matrix factorization and Krylov basis orthogonalization. We establish approximation error bounds for the Newton-L\'eja interpolation and derive an inexact oracle error bound that separates stochastic trace estimation from the approximation errors in the projected and residual components. Numerical experiments on both large sparse synthetic matrices (Gaussian Markov Random fields) and large-scale real-world matrices demonstrate that our method substantially reduces computational time while retaining competitive accuracy, thereby confirming its effectiveness, scalability, and favorable accuracy-cost trade-off for large-scale log-determinant estimation compared with stochastic Lanczos quadrature (SLQ) and related stochastic trace-estimation methods.