Complexity of Sparse Matrix Cholesky decomposition

This can only be answered exactly for abitrary matrices if P=NP ... so it's not possible to answer in general. The time complexity depends on the fill-reducing ordering used, which is attempting to get an approximate solution to an NP hard problem.

However, for the very special case of a matrix coming from a regular square 2D or 3D mesh, there is an answer. In this case, nested dissection gives an ordering that is asymptotically optimal. For a 2D s-by-s mesh, the matrix has dimension n = s^2 and I think about 5n entries. In this case, L has 31*(n log2(n)/8)+O(n) nonzeros, and the work is 829*(n^(3/2))/84+O(n log n). For a 3D s-by-s-by-s mesh with n = s^3, there are O(n^(4/3)) nonzeros in L and O(n^2) operations are required to compute L.

The Python Library Numpy (Numerical Python) also has a module for cholesky - np.linalg.cholesky, I have provided the link to the docs although I am not sure if this answers the question, might need some experimentation.

It's worth looking at an incomplete Cholesky decomposition, which there multiple variations of but typically either only compute the entries in the triangular factor that are nonzero in the input, or use a low rank approximation of the decomposition. It's not clear from your question why you are interested in the asymptotic complexity but with the Gaussian process tag I can guess that you're decomposing a covariance kernel matrix and repeatedly solving a linear system during inference. The incomplete factorization is most often used as a preconditioner in these kinds of applications- while it's not exact, it's very efficient, easy to update incrementally, and can greatly accelerate solvers.

In application, however, the answer to your question is very, very likely to be irrelevant to the approach that performs best for your purposes. There are $O(n log(n)^k)$ time algorithms for incomplete factorizations that are of no practical use for the same reasons as the algorithms for matrix multiplication with the lowest exponents. Knowing what precisely you need for your application will inform your options but you already have far and away the best tools to find out- taking a few minutes writing some code to generate synthetic data or sample varying sizes of your real data and timing it in the same way that you're interested in for application. If you're doing one factorization and solving many systems, the time to do the factorization is likely going to be dwarfed by the solving. If you are doing many factorizations, especially from scratch each time, the constant and linear factors of the running time will have even more impact. Plus, the sparsify patterns and kernel itself can have an enormous impact on the performance with the same algorithm. It's not difficult to construct kernels that have fully dense covariance matrices, tridiagonal precision matrices, and a Cholesky factorization that is just the lower triangular portion of the matrix exactly up to scaling by a diagonal matrix (an exponential kernel in 1d has all 3 simultaneously). Profile first, optimize last.