How for loops, map and comprehensions act on variables?

[Marco-Congedo]

I am writing an iterative algorithm. The function here below is part of the main function and performs a line search dividing by two the step size until the loss function is decreased. All variables are appropriately initialized in the main function and the algorithm works well. Now, if i switch line (1) with line (2) or (3), which are here below commented, the overall algorithm converges slower. Some debugging shows that the computation of the loss somehow changes and that this happens only if the for i loop does not break all the times at the first pass. I am puzzled. Any idea what is going here?

function _linesearch()
    for i = 1:lsmax
        M = (1.0/i * direction) + I
        B₊ = B * M
        𝐃₊ = [Hermitian(M'*D*M) for D ∈ 𝐃] # (1)
        # for j=1:length(𝐃) 𝐃₊[j] = Hermitian(M'*𝐃[j]*M) end # (2)
        # map!(D -> Hermitian(M'*D*M), 𝐃₊, 𝐃) # (3)
        loss₊ = -(logabsdet(B₊)[1]) + 0.5*sum(mean(log, [Diagonal(D) for D ∈ 𝐃₊])) 
        loss₊ < loss && break
    end
    return 𝐃₊, B₊, loss₊
end

PS: in the main function _linesearch is called as

𝐃, B, loss = _linesearch() 

and 𝐃 is used to find a new direction.

[ericphanson]

Is the only thing you change switching which of those lines are commented? Since versions (2) and (3) mutate 𝐃₊, which isn’t created within the body of the function _linesearch, meaning if you only comment (1) and uncomment e.g. (2), then they must be referring to a global variable 𝐃₊. That could explain some performance loss at least. (Version (1) on the other hand creates 𝐃₊ so it never refers to a global within that function).

[Marco-Congedo]

Yes, it is the only thing i do, just use one of the three options at the time.
I thought that if 𝐃₊ exists globally, then 𝐃₊ = [Hermitian(M'*D*M) for D ∈ 𝐃] writes into the existing variable. But in any case, even allowing this difference between line (1) versus either (2) or (3), i don’t get why in this code the computation of loss₊ gives different results with the different options.

To make things more uncomprehensible to me, if i use this line, as i said, it works well:

𝐃₊ = [Hermitian(M'*D*M) for D ∈ 𝐃]

if i use two lines one after the other like this

𝐃₊ = [Hermitian(M'*D*M) for D ∈ 𝐃]
for j=1:length(𝐃) 𝐃₊[j] = Hermitian(M'*𝐃[j]*M) end

it still works well and does the same thing. But if i use these two lines one after the other

for j=1:length(𝐃) 𝐃₊[j] = Hermitian(M'*𝐃[j]*M) end
𝐃₊ = [Hermitian(M'*D*M) for D ∈ 𝐃]

it gets skrewed up.
Why?

[stevengj]

@Marco-Congedo

All variables are appropriately initialized in the main function and the algorithm works well .

It looks like you are using global variables (your linesearch() function takes no arguments at all!)? Don’t.

This is bad for performance in Julia, and it is also bad for program structure (in any language) for a function to depend on (non-constant) global variables. Pass data into your function through parameters. If all inputs are passed to a function via parameters, and outputs are produced only via the return value (or mutating parameters), that makes the function self-contained (how it works doesn’t depend on the context in which is called) and hence much easier to re-use and understand.

If _linesearch() is an inner function referencing local variables from the surrounding scope, then you could be running into this problem: performance of captured variables in closures · Issue #15276 · JuliaLang/julia · GitHub, for which the workaround is to wrap captured variables in a let block. But it makes little sense to have an inner function that takes no parameters and is called only once — why have a function at all? (If by “main function” you mean that you wrapped your entire program in one giant function + inner functions, then you should re-think your code organization — this is hardly better than using global variables from a code-structure standpoint.)

[ericphanson]

So the line

𝐃₊ = [Hermitian(M'*D*M) for D ∈ 𝐃]

doesn’t mutate any preexisting variable 𝐃₊; instead, it points the “binding” 𝐃₊ at the vector [Hermitian(M'*D*M) for D ∈ 𝐃]. Inside a function you are in local scope, and a new local variable 𝐃₊ is created. To override that, you could write

global 𝐃₊
𝐃₊ = [Hermitian(M'*D*M) for D ∈ 𝐃]

which says that “in this function, I want 𝐃₊ to mean the global one, don’t create a new local one”.

The line

for j=1:length(𝐃) 𝐃₊[j] = Hermitian(M'*𝐃[j]*M) end

is necessarily mutating; 𝐃₊ must exist for 𝐃₊[j] to be accessible, and 𝐃₊[j] = Hermitian(M'*𝐃[j]*M) says “please change the existing thing 𝐃₊[j]”. Since there is no local 𝐃₊, the only way for that work is if it modifies the global 𝐃₊, so it does so. See the manual section on scoping for more. You can test this with a simpler example:

Case 1: modify the global

julia> D = [0]
1-element Array{Int64,1}:
 0

julia> function test0()
       for j = eachindex(D)
              D[j] = j
       end
       return D
       end
test0 (generic function with 1 method)

julia> test0()
1-element Array{Int64,1}:
 1

julia> D
1-element Array{Int64,1}:
 1

Case 2: create a new local binding shadowing the global:

julia> D = [0]
1-element Array{Int64,1}:
 0

julia> function test1()
            D = [5]
            for j = eachindex(D)
                D[j] = j
           end
          return D
       end
test1 (generic function with 1 method)

julia> test1()
1-element Array{Int64,1}:
 1

julia> D
1-element Array{Int64,1}:
 0

For your second question, in the first version you create a local 𝐃₊, then mutate it. In the second version, that should just error because the line 𝐃₊ = ... should create a local binding (no matter in what order it’s coded in the function, i.e. before or after the for loop), and then the for loop line is then trying to modify the local binding before it’s assigned. E.g.

julia> function test2()
           for j = eachindex(D)
               D[j] = j
           end
           D = [5]
           return D
       end
test2 (generic function with 1 method)

julia> test2()
ERROR: UndefVarError: D not defined
Stacktrace:
 [1] test2() at ./REPL[19]:2
 [2] top-level scope at REPL[20]:1 

edit: I very much agree with @stevengj’s post above; the solution here is to not use global variables. I just wanted to take the opportunity to explain what’s going on, because the scoping can be a bit subtle but once you understand it, it all makes sense.

[Marco-Congedo]

thank tou @ stevengj for the suggestion and @ ericphanson for the lengthy explanations. It starts making sense to me. I am a confused about the suggestion “not use global variables”.

The _linesearch is called by the main function at each iteration and returns the updated values that are needed. I though that it is a good practice to create memory for all variables that are used over and over again at the beginning of the function and then allowing the function to manipulate those, reusing the same allocated memory, as compared to create again and again variables on the fly as they are needed. If so, each time _linesearch is called i would prefer to write into an existing global 𝐃₊ vector (i.e., using the for loop line in my code), rather then creating one (i.e., the comprehension generator line in my code). Is that correct?

[stevengj]

@Marco-Congedo

I though that it is a good practice to create memory for all variables that are used over and over again at the beginning of the function and then allowing the function to manipulate those, reusing the same allocated memory, as compared to create again and again variables on the fly as they are needed.

Preallocating arrays is often a good practice, but you can still pass them as arguments — it doesn’t require you to use globals.

In any case, as @ericphanson pointed out above, a statement like B₊ = B * M or 𝐃₊ = [Hermitian(M'*D*M) for D ∈ 𝐃] allocates new arrays anyway — you are “rebinding” the variables B₊ and 𝐃₊ to “point” to new arrays in memory, not mutating them in-place. See also this discussion.

(This is not unique to Julia! If you do a = [3,4,5], then b = a, followed by a = [4,5,6] in Python or Matlab or most other languages with analogous operations, b will still be [3,4,5].)

[Marco-Congedo]

Ok, thanks a lot.

[DNF]

It seems like you have a bit of a wrong perspective on how bindings work. This is an old, but good explanation of the sort of issues you seem to have: Values vs. Bindings: The Map is Not the Territory · John Myles White

The explanation uses Julia code which may in some cases be a bit out of date.

[kevbonham]

I can’t add much of substance, other than a note that in julia, the convention is that functions that are mutating usually have ! in the function signature. So if you end up re-writing _linesearch() to take an argument and then mutate it the convention would be to name it _linesearch!(D) or whatever.

Knowing this can also be a hint about which functions will operate in-place vs allocate a new vector. For example, compare

julia> x = Vector(1:4);

julia> replace(x) do i
           isodd(i) ? i : 0
       end
4-element Array{Int64,1}:
 1
 0
 3
 0

julia> x
4-element Array{Int64,1}:
 1
 2
 3
 4

to

julia> replace!(x) do i
           isodd(i) ? i : 0
       end
4-element Array{Int64,1}:
 1
 0
 3
 0

julia> x
4-element Array{Int64,1}:
 1
 0
 3
 0

[Marco-Congedo]

@ DNF thanks, very useful indeed. I think i found was i was looking for: a syntax to write into an existing matrix (like A[:]=B*C), not to create a new matrix (like A=B*C). By the way, what is the difference between A[:]=B*C and A.=B*C ?