# Very long time for addition of complex identity matrix to transposed sparse matrix

**URL:** https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465
**Category:** Performance
**Tags:** question, sparsearrays
**Created:** [May 30, 2025, 10:18am UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465 "2025-05-30T10:18:19Z")
**Posts on this page:** 11
**Page:** 1

<div class="post-metadata">

### Author: ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 30, 2025, 10:18am UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/1 "2025-05-30T10:18:19Z")

</div>

Hi,

I am surprised by these timings on 1.10. Can it be improved?

```julia
using SparseArrays, LinearAlgebra
A = sprandn(90000,90000, 1e-5)
@time (A + complex(0,1.) * I); # 0.001344 seconds (13 allocations: 7.342 MiB)
@time (transpose(A) + complex(0,1.) * I); #101.973608 seconds (193.74 k allocations: 21.416 MiB, 0.08% compilation time)

```

---

<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: [May 30, 2025, 12:41pm UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/2 "2025-05-30T12:41:07Z")

</div>

`A + x*I` is hitting a [specialized method](https://github.com/JuliaSparse/SparseArrays.jl/blob/6d072a81fca5f4394f88a012f4ce914c70769303/src/sparsematrix.jl#L4388-L4395), whereas `A' + x*I` hits a [generic method](https://github.com/JuliaLang/LinearAlgebra.jl/blob/c4a477fff61b1f72c7d867e3eebad8cc1eab4d8b/src/uniformscaling.jl#L220-L227) that mutates one element of the diagonal at a time, which is very slow for sparse matrices.

Should be easily fixable by adding some specialized methods, e.g.

```julia
import Base: +, -
using LinearAlgebra: AdjOrTrans, wrapperop
(+)(A::AdjOrTrans{<:Any, <:AbstractSparseMatrix}, J::UniformScaling{<:Number}) =
    wrapperop(A)(parent(A) + wrapperop(A)(J))
(-)(A::AdjOrTrans{<:Any, <:AbstractSparseMatrix}, J::UniformScaling{<:Number}) =
    wrapperop(A)(parent(A) - wrapperop(A)(J))
(+)(J::UniformScaling{<:Number}, A::AdjOrTrans{<:Any, <:AbstractSparseMatrix}) =
    wrapperop(A)(wrapperop(A)(J) + parent(A))
(-)(J::UniformScaling{<:Number}, A::AdjOrTrans{<:Any, <:AbstractSparseMatrix}) =
    wrapperop(A)(wrapperop(A)(J) - parent(A))

```

Should be an easy PR for someone?

---

<div class="post-metadata">

### Author: ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)
#### Post date: [May 30, 2025, 12:43pm UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/3 "2025-05-30T12:43:14Z")

</div>

> <https://github.com/JuliaLang/LinearAlgebra.jl/pull/1366>
>
> Avoids hitting slow generic methods.
> Discussed in https://discourse.julialang.o…rg/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465.
> After this, the following are comparable:
> \`\`\`julia
> julia\> A = sprandn(90, 90, 0.01);
> 
> julia\> @btime ($A + complex(0,1.) \* I);
> 1.313 μs (13 allocations: 8.02 KiB)
> 
> julia\> @btime (transpose($A) + complex(0,1.) \* I);
> 1.416 μs (13 allocations: 8.02 KiB)
> \`\`\`

---

<div class="post-metadata">

### Author: ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 31, 2025, 7:32am UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/4 "2025-05-31T07:32:16Z")

</div>

Is it a similar issue for this one?

```julia
A = sprandn(9000,9000, 1e-5)
@time (A + complex(0,1.) * I); # 0.000110 seconds (13 allocations: 581.859 KiB)
A = @view sprandn(9000,9000, 1e-5)[1:end-1, 1:end-1]
@time (A + complex(0,1.) * I); # 0.351283 seconds (13 allocations: 539.094 KiB)

```

---

<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: [May 31, 2025, 7:47am UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/5 "2025-05-31T07:47:42Z")

</div>

I think this one is just a version of “views on sparse arrays suck”

---

<div class="post-metadata">

### Author: ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 31, 2025, 7:59am UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/6 "2025-05-31T07:59:45Z")

</div>

It took me a long time to understand why my bifurcation code for some PDE was really slow ☹  
This tread shows why. Do you know of a package which solves this?

---

<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: [May 31, 2025, 8:01am UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/7 "2025-05-31T08:01:45Z")

</div>

I think the problem with the current sparse view is that it creates a wrapper where you can access individual matrix elements, but these accesses are very slow for sparse matrices (each one requires a `searchsorted`), and other operations beyond `getindex(A, i, j)` are not well-optimized.

As a result, whenever I need a view of a sparse matrix, I build it myself using the fields of the `SparseMatrixCSC` type. Once you get the hang of it, it’s not overly difficult.

---

<div class="post-metadata">

### Author: ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 31, 2025, 8:13am UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/8 "2025-05-31T08:13:26Z")

</div>

In `BifurcationKit`, I have in many places `@view A ....` where A is an AbstractMatrix.  
Basically, I do a lot of `sigma * I + B` where B can be a view to a sparse matrix. Maybe I should code this differently since it allocates anyway

---

<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: [May 31, 2025, 2:52pm UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/9 "2025-05-31T14:52:20Z")

</div>

Can you explain a little bit more about the context in which you need these views?

---

<div class="post-metadata">

### Author: ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 31, 2025, 4:56pm UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/10 "2025-05-31T16:56:10Z")

</div>

It happened during fold continuation as in [here](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/codim2Continuation/#Fold-continuation-(theory)). I need to form the matrix M\_f as jacobian for newton iterations. Later, I need to compute the spectrum of `dF` so I form `dF = @view Mf[1:end-1,1:end-1]` and compute the spectrum. Unfortunately, `dF` is not well conditioned, so I use a Shift-inverse strategy which requires to precompute `factorize(sigma * I - dF)`. In a nutshell, all the problems of this post appears.

> I know I dont need to form the bordered system to solve Mf, I could use a bordering strategy, it just happens to be slower in my use case.

If I were doing Hopf continuation, it would be the same. Actually, for PD/NS continuation (see [here](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/PDContinuationPO/#Algorithmic-details)) it would be similar.

---

<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: [May 31, 2025, 5:20pm UTC](https://discourse.julialang.org/t/very-long-time-for-addition-of-complex-identity-matrix-to-transposed-sparse-matrix/129465/11 "2025-05-31T17:20:12Z")

</div>

Maybe you could look at it the other way, and consider the `Mf[1:end-1,1:end-1]` matrix to be your starting point, then either create a bigger block matrix or perform a manual solve for the system of the Newton iteration? Idk if there are formulas for solving \begin{pmatrix} M & u \\ v^T & 0 \end{pmatrix} = \begin{pmatrix} a \\ b \end{pmatrix} directly?

Alternately you could keep both `M` and its principal submatrix in storage, and perform the update manually with a dedicated version in case `M` is sparse?
