Running out of RAM with repeated cholesky solves

Specific Domains Numerics
,
1

Hello,

I have a system in which I need to repeatedly solve an electrostatics problem with an evolving charge configuration. When I left my program running overnight, I came back to find that RAM utilization had accumulated to around 30GB and the program had slowed to a crawl. I believe I’ve located the problem line
ldiv!(solution, stiff_factorization, force_vector)
When I comment this out, the program is able to complete with well under ~1GB of RAM allocation. I’m wondering if I need to manually clean up any memory being used by the solver. I tried to create a small example of what I’m working with, but memory doesn’t seem to accumulate in the same way. I’m not very familiar with the workings of the Garbage Collector, so any advice about how to avoid a memory leak over longer runs is appreciated. Thank you!

using JLD2
using LinearAlgebra

stiffness_matrix = load("./dev/ldiv_memory/stiff_mat.jld2", "stiff_mat")
issymmetric(stiffness_matrix)
factorization = cholesky(stiffness_matrix)

dim = size(stiffness_matrix, 1)
sol = zeros(dim)

for i=1:1_000_000
    force_vector = rand(dim)
    ldiv!(sol, factorization, force_vector)
end
2

What Julia version are you using. I fixed something similar and it landed on 1.12.5 so if you could try there

1 Like
3

I just ran juliaup update and up in the package manager. I’m on the 1.12.5 release, but task manager (Windows) keeps showing the memory increasing. Previously, I tried using cg! which seemed to reduce the memory load but it also seems to be running more slowly for my matrix.

4

You really shouldn’t need to manually clean up memory. In the Windows task manager, what memory metric are you looking at? (Depending on how you measure memory it can be deceptive due to virtual memory.)

Is stiffness_matrix a dense matrix or a sparse one?

(In the latter case, it might be a memory leak in the external CHOLMOD library that we are calling for sparse Cholesky solves? Though this would be surprising… that library has been around for a long while.)

Are you sure it’s not something else? For example, if NaN’s have crept into your vectors then things will slow to a crawl due to floating-point exceptions.

Does the example code you posted exhibit the problem (if you use it with your stiff_mat.jld2) or not?

1 Like
5

Thanks for the feedback, it led me down some interesting investigations. My stiffness matrix is assembled as a dimension ~5000 sparse matrix SparseMatrixCSC{Float64, Int64} but I had cast it to a dense matrix when I saved it. I’ve modified the snippet to cast to a sparse matrix, and now I see a steady memory usage increase that drops down after I terminated the program after 665800 iterations.

image
image799×1013 20 KB

My full program drastically slows down after the memory allocation reaches 80%, which is why I suspect that it’s not just the way that task manager measures memory usage. I’m using the solution vector for downstream calculations, which haven’t thrown an errors for values being too large (I can try to re-run with an explicit NaN check as well).

I tried to create a version of the code that doesn’t need the file to run (Also happy to provide the saved matrix but it seems I’m not allowed to upload it here). The synthetic random matrix takes a lot longer to churn through iterations but I think I see memory uptick as well

using JLD2
using LinearAlgebra
using SparseArrays

load_file = false
if load_file
    stiffness_matrix = load("./stiff_mat.jld2", "stiff_mat")
    println(typeof(stiffness_matrix))
    stiffness_matrix = SparseMatrixCSC(stiffness_matrix)
    println(typeof(stiffness_matrix))
else
    mock_size = 5000
    mock_density = 0.002
    A = sprand(Float64, mock_size, mock_size, mock_density)
    S = A + A'
    
    # 3. Make it strictly diagonally dominant to guarantee Positive Definiteness
    # Since S has positive entries, sum() gives the sum of absolute off-diagonals
    row_sums = vec(sum(S, dims=2))
    
    # Create a diagonal matrix where each diagonal is slightly larger than the row sum
    D = spdiagm(0 => row_sums .+ 1.0)
    
    stiffness_matrix = SparseMatrixCSC{Float64, Int64}(S + D)
end
println("Is matrix symmetric: $(issymmetric(stiffness_matrix)), determinant: $(det(stiffness_matrix))")
factorization = cholesky(stiffness_matrix)

dim = size(stiffness_matrix, 1)
sol = zeros(dim)

for i=1:1_000_000
    force_vector = rand(dim)
    if i%100 == 0
        println("Iteration $(i)")
    end
    ldiv!(sol, factorization, force_vector)
    if !all(isfinite.(sol))
        throw(ErrorException("Got a non-finite solution"))
    end
end
6

You can upload it elsewhere (e.g. google drive, github gist, dropbox, …) and post a link here.

7

Uploaded to Google drive: https://drive.google.com/file/d/1GyU1ZflSY9-WyAMQKhVQu0me2mUIK_IO/view?usp=drive_link