# Best way to repeatedly multiply a Vector by a SparseMatrix

**URL:** <https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355>\
**Category:** General Usage\
**Created:** [August 31, 2018, 8:08am UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355 "2018-08-31T08:08:19Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![teored90](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/teored90/32/5142_2.png) [@teored90](https://discourse.julialang.org/u/teored90)\
**Post date:** [August 31, 2018, 8:08am UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/1 "2018-08-31T08:08:19Z")

</div>

Hi, in my code I need to evaluate A^n \* v where A is a possibly very large sparse matrix and v is a vector.

I suppose it is much faster to perform `A * (A * ... * (A * v)` rather than elevating A to the n-th power, but what is the most efficient way to do so?

I currently have

```julia
function A_pow_n_B!(res, A, n::Integer, v)
    tmp = similar(v)
    mul!(res, A, v)
    for i in 1:n-1
        mul!(tmp, A, res)
        copyto!(res, tmp)
    end
end

```

Can this code be further improved?

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [August 31, 2018, 10:04am UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/2 "2018-08-31T10:04:33Z")

</div>

I guess so, as the `copyto!` isn’t necessary in every iteration. You could alternate between `mul!(tmp, A, res)` and `mul!(res, A, tmp)` and return `tmp` or `res` in the end depending on which holds the result.

---

<div class="post-metadata">

**Author:** ![teored90](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/teored90/32/5142_2.png) [@teored90](https://discourse.julialang.org/u/teored90)\
**Post date:** [August 31, 2018, 11:52am UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/3 "2018-08-31T11:52:28Z")

</div>

What you say makes a lot of sense.  
But I tried to do as you suggested but the time difference seems negligible (most of the time is spent doing `mul!`: in my use case n \sim 10)

```julia
function new_A_pow_n_B!(res, x, n::Integer, y)
    tmp = similar(y)
    if isodd(n)
        mul!(res, x, y)
        for i in 2:2:n
            mul!(tmp, x, res)
            mul!(res, x, tmp)
        end
    else
        mul!(tmp, x, y)
        for i in 2:2:n-1
            mul!(res, x, tmp)
            mul!(tmp, x, res)
        end
        mul!(res, x, tmp)
    end
end

```

I get

```julia
d = 5000
r = zeros(d)

@benchmark new_A_pow_n_B!($(r), $(sprand(ComplexF64, d, d, .01)), 9, $(rand(ComplexF64, d)))
BenchmarkTools.Trial: 
  memory estimate: 78.20 KiB
  allocs estimate: 2
  --------------
  minimum time: 3.181 ms (0.00% GC)
  median time: 3.551 ms (0.00% GC)
  mean time: 3.696 ms (0.05% GC)
  maximum time: 7.348 ms (0.00% GC)
  --------------
  samples: 1349
  evals/sample: 1

@benchmark A_pow_n_B!($(r), $(sprand(ComplexF64, d, d, .01)), 9, $(rand(ComplexF64, d)))
BenchmarkTools.Trial: 
  memory estimate: 78.20 KiB
  allocs estimate: 2
  --------------
  minimum time: 3.254 ms (0.00% GC)
  median time: 3.716 ms (0.00% GC)
  mean time: 3.947 ms (1.17% GC)
  maximum time: 53.210 ms (87.77% GC)
  --------------
  samples: 1264
  evals/sample: 1

```

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [August 31, 2018, 11:58am UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/4 "2018-08-31T11:58:07Z")

</div>

FWIW, the gain should depend on `d` though. So for very large `d` … 🙂

---

<div class="post-metadata">

**Author:** ![RaulDurand](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rauldurand/32/2044_2.png) [@RaulDurand](https://discourse.julialang.org/u/RaulDurand)\
**Post date:** [August 31, 2018, 12:09pm UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/5 "2018-08-31T12:09:53Z")

</div>

Singular value decomposition can be used to calculate A^n in one step. Depending on the size of A and n it may be cheaper.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [August 31, 2018, 12:19pm UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/6 "2018-08-31T12:19:54Z")

</div>

> [@RaulDurand](#):
>
> Depending on the size of A and n it may be cheaper.

… and have smaller error, also depending on `A` and `n`.

---

<div class="post-metadata">

**Author:** ![teored90](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/teored90/32/5142_2.png) [@teored90](https://discourse.julialang.org/u/teored90)\
**Post date:** [August 31, 2018, 12:21pm UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/7 "2018-08-31T12:21:00Z")

</div>

There is no requirement on the size, but let’s say it can vary from 100 to 10^6, while n \lesssim 10. The matrix is sparse. Do you think it is worth to consider SVD?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [August 31, 2018, 12:47pm UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/8 "2018-08-31T12:47:41Z")

</div>

> [@teored90](#):
>
> The matrix is sparse.

What people call “sparse” is very subjective, so it depends. But perhaps for elements mostly around the diagonal and no tricky conditioning/cancellation issues, I would just go for the algorithm you first posted; as SVD may be unnecessarily complex. There are partial algorithms like

> **[GitHub - bwlewis/irlba: Fast truncated singular value decompositions](https://github.com/bwlewis/irlba)**
>
> Fast truncated singular value decompositions. Contribute to bwlewis/irlba development by creating an account on GitHub.

but this may be overkill.

---

<div class="post-metadata">

**Author:** ![teored90](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/teored90/32/5142_2.png) [@teored90](https://discourse.julialang.org/u/teored90)\
**Post date:** [August 31, 2018, 3:44pm UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/9 "2018-08-31T15:44:54Z")

</div>

The thing is that I can’t assume anything about the structure of the matrix, or it’s density. But, sticking to the repeated multiplication in the OP, are there other optimizations possible with Julia?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [August 31, 2018, 4:09pm UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/10 "2018-08-31T16:09:32Z")

</div>

I think your original function is fine, I would just make sure that `v` is not sparse to begin with if I don’t expect sparsity to be preserved. But if you literally can’t assume anything, then it is hard to say more.

You can also optimize multiplication (since that’s the only thing that matters). Transposing `A` to `B` (in a wrapper function) and then multiplying by `B'` may speed up things a bit, because of contiguity.

---

<div class="post-metadata">

**Author:** ![Nectarineimp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nectarineimp/32/883_2.png) [@Nectarineimp](https://discourse.julialang.org/u/Nectarineimp)\
**Post date:** [August 31, 2018, 9:03pm UTC](https://discourse.julialang.org/t/best-way-to-repeatedly-multiply-a-vector-by-a-sparsematrix/14355/11 "2018-08-31T21:03:35Z")

</div>

Maybe this is dumb, but instead of iterating the entire giant matrix, why not just get the coordinates of the structurally non-zero elements and operate on just those?

[`SparseArrays.findnz`](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#SparseArrays.findnz)