# Help on optimization

**URL:** <https://discourse.julialang.org/t/help-on-optimization/10238>\
**Category:** General Usage\
**Created:** [April 8, 2018, 10:10pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238 "2018-04-08T22:10:34Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [April 8, 2018, 10:10pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/1 "2018-04-08T22:10:34Z")

</div>

Hello everyone,

I need some help on optimization of a function, which looks like this

```julia
# re-edit
SpMatrixListUnion = Union{Vector{SparseMatrixCSC{Float64,Int64}},Vector{SparseMatrixCSC{Complex128,Int64}}}

function avg_V_on_ψ!(
    res::Array{Float64,2}, # result to write
    vec_list::Vector{Matrix{Float64}}, # a list of V's, 
    OP_list::SpMatrixListUnion # a list of O's
    ) # the computation is V' * O * V, 
                                         #result is written into Array "res" in a consistent manner
    vl = length(vec_list)
    opl = length(OP_list)
    for i=1:vl
        I,J,V = findnz(OP_list[i]) # helpful when OP_list[i] is a very sparse matrix
        res[:,i] = real.(sum((vec_list[i][I,:]).*broadcast(*,V,vec_list[i][J,:]),1)) # V' * O * V for sparse matrix O
    end
    return nothing
end

```

For each matrix `O=OP_list[i]` in OP\_list and each vector set `V=vec_list[i] in vec_list`, compute the average `V'*O*V` and write result to res. I am working with sparse matrices for O in OP\_list, so I have done a first step optimization by using findnz. For really sparse matrices, the result is good but not satisfactory.

It seems to me that there cannot be further reduction on the number of arithmatic operations. So I hope to ask some experienced users for help, to make such function even faster, because it is the only bottleneck in my program …

Thanks in advance in any suggestions!

BTW you can verify the code by running the following test:

```julia
# avg_V_on_ψ0 is a slow but guaranteed version for the task

avg_V_on_ψ0(vecs,OP) = [real(dot(vecs[:,k],OP*vecs[:,k])) for k=1:size(vecs)[2]]

# the test

vlist = [rand(400,400) for k=1:300];
oplist = [1000.0sprand(400,400,0.0001) for k=1:300];
res2 = zeros(Float64,400,300);
gc()
@time res1 = [avg_V_on_ψ0(vlist[i],oplist[i]) for i=1:300];
@time avg_V_on_ψ!(res2,vlist,oplist)

# test : norm should be equal ~
norm(hcat(res1...)-res2)

```

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [April 9, 2018, 12:26am UTC](https://discourse.julialang.org/t/help-on-optimization/10238/2 "2018-04-09T00:26:22Z")

</div>

Would you mind [quoting your code](https://discourse.julialang.org/t/psa-how-to-quote-code-with-backticks/7530) to make it more readable?

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [April 9, 2018, 1:48am UTC](https://discourse.julialang.org/t/help-on-optimization/10238/3 "2018-04-09T01:48:03Z")

</div>

Your code doesn’t run here. I’m assuming that `SpMatrixUnion` should be `SpMatrixListUnion`, but then I still get

```julia
julia> avg_V_on_ψ!(res2,vlist,oplist)
ERROR: MethodError: no method matching avg_V_on_ψ!(::Array{Float64,2}, ::Array{Array{Float64,2},1}, ::Array{SparseMatrixCSC{Float64,Int64},1})
Closest candidates are:
  avg_V_on_ψ!(::Array{Float64,2}, ::Array{Array{Float64,2},1}, ::Union{SparseMatrixCSC{Complex{Float64},Int64}, SparseMatrixCSC{Float64,Int64}}) at REPL[4]:6

```

on 0.6.2.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 9, 2018, 8:54am UTC](https://discourse.julialang.org/t/help-on-optimization/10238/4 "2018-04-09T08:54:40Z")

</div>

Using views give a meagre 20% improvement:

```julia
function avg_V_on_ψ!2(
    res::Array{Float64,2},
    vec_list::Vector{Matrix{Float64}},
    OP_list
    )
    for i=1:length(vec_list)
        I,J,V = findnz(OP_list[i])
        v = vec_list[i]
        vI = view(v, I, :)
        vJ = view(v, J, :)
        res[:,i] = real.(sum(vI .* broadcast(*, V, vJ), 1))
    end
    return nothing
end

```

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [April 9, 2018, 2:00pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/5 "2018-04-09T14:00:52Z")

</div>

Sorry I’ve got a paranoid way of managing types. I have replaced the `SpMatrixUnion` by `SpMatrixListUnion`. This time it should work …

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [April 9, 2018, 2:04pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/6 "2018-04-09T14:04:00Z")

</div>

Excellent ! Thanks a lot.  
It seems that the `broadcast` function prefers views of matrices, instead of matrices themselves.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 9, 2018, 2:06pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/7 "2018-04-09T14:06:46Z")

</div>

I think the performance gain is just by avoiding having to copy the sub-matrices.

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [April 9, 2018, 2:14pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/8 "2018-04-09T14:14:15Z")

</div>

BTW kristofer.carlsson just removed the type after `OP_list`, that makes the code more transparent …

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [April 9, 2018, 2:27pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/9 "2018-04-09T14:27:53Z")

</div>

I just realized that `real.(sum(vI .* broadcast(*, V, vJ), 1))` is equivalent to `res[:,i] = real.(RowVector(V)*(vI.*vJ))`. The latter leads to 20% gain

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 9, 2018, 2:31pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/10 "2018-04-09T14:31:24Z")

</div>

Or `V' * (vI .* vJ)` (a bit more compact)

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [April 9, 2018, 2:34pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/11 "2018-04-09T14:34:33Z")

</div>

> [@kristoffer.carlsson](#):
>
> for i=1:length(vec\_list)  
> I,J,V = findnz(OP\_list[i])  
> v = vec\_list[i]  
> vI = view(v, I, 🙂  
> vJ = view(v, J, 🙂  
> res[:,i] = real.(sum(vI .\* broadcast(\*, V, vJ), 1))  
> end

yes, only if the O in op\_list has all real entries

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [April 9, 2018, 2:39pm UTC](https://discourse.julialang.org/t/help-on-optimization/10238/12 "2018-04-09T14:39:49Z")

</div>

> [@Lian\_Yunlong](#):
>
> I,J,V = findnz(OP\_list[i]) # helpful when OP\_list[i] is a very sparse matrix  
> res[:,i] = real.(sum((vec\_list[i][I,:])._broadcast(_,V,vec\_list[i][J,:]),1)) # V’ \* O \* V for sparse matrix O  
> end

optimized code

```julia
function avg_V_on_ψ!(
    res::Array{Float64,2},
    vec_list::Vector{Matrix{Float64}},
    OP_list::SpMatrixListUnion
    )
    """ for each operator in OP_list and each eigenvector set in vec_list, compute the average and write result to res."""
    for i=1:length(vec_list)
        I,J,V = findnz(OP_list[i])
        v = vec_list[i]
        vI = view(v, I, :)
        vJ = view(v, J, :)
        res[:,i] = real.(RowVector(V)*(vI.*vJ))
    end
    return nothing
end

```

before optimzation  
0.480163 seconds (59.48 k allocations: 607.347 MiB, 10.69% gc time)  
after  
0.331370 seconds (50.49 k allocations: 177.589 MiB, 4.70% gc time)