# How to efficiently update a Sparse Matrix (Diagonal + Sparse) in-place inside a hot loop?

**URL:** <https://discourse.julialang.org/t/how-to-efficiently-update-a-sparse-matrix-diagonal-sparse-in-place-inside-a-hot-loop/134705>\
**Category:** General Usage\
**Tags:** question, memory-allocation, sparse\
**Created:** [December 23, 2025, 10:10am UTC](https://discourse.julialang.org/t/how-to-efficiently-update-a-sparse-matrix-diagonal-sparse-in-place-inside-a-hot-loop/134705 "2025-12-23T10:10:46Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Matt\_jl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matt_jl/32/52364_2.png) [@Matt\_jl](https://discourse.julialang.org/u/Matt_jl)\
**Post date:** [December 23, 2025, 10:10am UTC](https://discourse.julialang.org/t/how-to-efficiently-update-a-sparse-matrix-diagonal-sparse-in-place-inside-a-hot-loop/134705/1 "2025-12-23T10:10:46Z")

</div>

Hi everyone,

I am working on an iterative algorithm where I need to construct and solve a linear system inside a “hot loop”.

Currently, my code allocates a new sparse matrix `A` at every iteration because of the operation between a `Diagonal` matrix and a `SparseMatrix`.

something like:

```julia-auto
using LinearAlgebra, SparseArrays, BenchmarkTools

function mwe(y::Vector{T}, max_iter::Int=10) where T
    m = length(y)
    
    D = spdiagm(0 => -2.0 * ones(T, m), 1 => ones(T, m-1), -1 => ones(T, m-1))
    AA = D' * D
    
    # Buffers
    w = ones(T, m)
    tmp_rhs = similar(y)
    result = similar(y)
    residuals = similar(y)

    for _ in 1:max_iter
        
        A = Diagonal(w) * AA

       
        tmp_rhs .= w .* y
        result .= A \ tmp_rhs

       # not used in mwe
        residuals .= y .- result
        w .= 1.0 ./ (1.0 .+ abs.(residuals)) 
    end
    
    return result
end

```

Bonus point: is possible to resolve the linear sistem in-place into the results vector ?  
How can I do that?

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [December 23, 2025, 11:09am UTC](https://discourse.julialang.org/t/how-to-efficiently-update-a-sparse-matrix-diagonal-sparse-in-place-inside-a-hot-loop/134705/2 "2025-12-23T11:09:25Z")

</div>

Just to clarify, inside your code I see a multiplication between diagonal and sparse, not an addition. Which is it?

---

<div class="post-metadata">

**Author:** ![Matt\_jl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matt_jl/32/52364_2.png) [@Matt\_jl](https://discourse.julialang.org/u/Matt_jl)\
**Post date:** [December 23, 2025, 2:24pm UTC](https://discourse.julialang.org/t/how-to-efficiently-update-a-sparse-matrix-diagonal-sparse-in-place-inside-a-hot-loop/134705/3 "2025-12-23T14:24:33Z")

</div>

typo sorry, I mean the multiplication.

But, mainly, I would like to solve the linear sistem in place directly into results vector.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 23, 2025, 3:51pm UTC](https://discourse.julialang.org/t/how-to-efficiently-update-a-sparse-matrix-diagonal-sparse-in-place-inside-a-hot-loop/134705/4 "2025-12-23T15:51:57Z")

</div>

You can multiply by a diagonal matrix in-place by looping over `rowvals` and `nonzeros`, just multiplying each element by the diagonal element of the corresponding row. In particular, allocate `A = copy(AA)` outside of the loop, and then do:

```julia-auto
for i in eachindex(nonzeros(A))
    nonzeros(A)[i] = w[rowvals(A)[i]] * nonzeros(AA)[i]
end

```

or, more compactly:

```julia-auto
nonzeros(A) .= getindex.(Ref(w), rowvals(A)) .* nonzeros(AA) # A = Diagonal(w) * AA

```

> [@Matt\_jl](#):
>
> Bonus point: is possible to resolve the linear sistem in-place into the results vector

You need to:

1. Allocate the LU factorization `F = lu(A)` outside of the loop
2. Re-use this to compute `lu!(F, A)` in-place in the loop (still has a few allocations)
3. use `ldiv!(result, F, tmp_rhs)` to use the factorization to solve the system in-place.

However, if you care about efficiency, there are probably other things that will have a bigger impact. Allocations aren’t usually the slowest part of solving large sparse systems. Instead, I would think about the following issues:

1. Your matrix A is banded. You might want to exploit this structure instead of using generic sparse matrices, by using BandedMatrices.jl
2. You are computing `result = (W * AA) \ (W * y)` where `W = Diagonal(w)`. But this is equivalent to computing `AA \ y`. Why are you multiplying by `W` at all? Is this a typo?
3. Assuming (2) is a typo, you still don’t need to update your sparse matrix. (W \texttt{AA})^{-1} = \texttt{AA}^{-1} W^{-1}, so you can simply divide by W first (equivalent to dividing elementwise by `w`) and then left-divide by `AA`. Since `AA` is not changing at all in your loop, you can compute its LU factorization `F = lu(AA)` _once_ and re-use it.
4. The matrix `AA = D' * D` is symmetric positive-definite, so rather than LU factorization, you can take advantage of this structure by `AA = Symmetric(D' * D)` and then using the `F = cholesky(AA)` factorization instead (more than twice as fast).

When you are doing large-scale linear algebra, you mainly want to think about the overall algebraic structure and the linear algebra, rather than trying to micro-optimize your loops and allocations.