In the spirit of @brrrt
(see #random > Rename LoopVectorization.@avx
?), I recently spent days debugging the following bug:
# get hidden neuron firing rate
reservoir.ξhidden!(cache.r)
@avx cache.r .+= reservoir.f.(u)
# get output neuron firing rate
reservoir.ξ!(cache.z)
if isexploring[]
reservoir.ξ!(cache.z)
@avx cache.z .+= reservoir.Wout * cache.r
else
@avx cache.z .= reservoir.Wout * cache.r
end
# update du
@avx du .= (-u .+ reservoir.λ * reservoir.Wr * cache.r .+
reservoir.Win * input(t) .+
reservoir.Wfb * cache.z) ./ reservoir.τ
This code fails to produce anything resembling working. A red herring related to the DE stepping is why I spent so many days on it until I found out that removing @avx
from the line corresponding to du
fixed things. Now, I am wondering why the error occurred (like a more detailed explanation than LoopVectorization.jl is unsafe). The size of the weight matrices are:
size(Wr) == (1000, 1000)
size(Win) == (1000, 1)
size(Wfb) == (1000, 1)
size(Wout) == (1, 1000)
(notice that Wout
is not part of a broadcast, it is a matmul)I think it's probably this:
image.png
Dimensions of size 1 are not kosher with LV, but I never really understood why
Yeah but that warning only mentions broadcasting. The full expression is z .= z .+ Wout * r
, so the 1x1000 array is embedded in an expression that uses broadcast, but it is not part of the broadcast itself. This is partially why I asked the question — to get clarity on if this kind of expression is still bad. And if the following would "fix" things:
tmp = Wout * r
@avx z .+= tmp
Also it is weird that I only needed to remove @avx
from the du
expression which involves no 1xN arrays!
What's the size of r
?
If r
is a 1000xN
array, then Wout * r
is a 1xN
array so the problem is not solved
You can make Wout
a Transpose{Vector}
though
Ah, you're right! r
is 1000x1, but I assume 1x1 arrays also not supported
Yeah let me try with the transpose
Thanks, this makes sense now
Is r
a Matrix
, or is it a Vector
?
You'll need r
to be a Vector
in order to not get a 1x1
array out
julia> let Wout = rand(1000)', r = rand(1000, 1)
Wout * r
end
1×1 adjoint(::Vector{Float64}) with eltype Float64:
244.01241695906867
julia> let Wout = rand(1000)', r = rand(1000)
Wout * r
end
236.71936258108883
It's a Vector
Last updated: Nov 22 2024 at 04:41 UTC