Assigning to a view of an array

New to Julia
,
1

The purpose of creating a view from an array is to avoid a copy.

However, a view cannot be assigned to.

julia> X = rand(3,3)
3×3 Matrix{Float64}:
0.135318  0.958828  0.0860895
0.715183  0.136236  0.542946
0.513489  0.974295  0.155925

julia> @views sX = X[1:2,1:2]
2×2 view(::Matrix{Float64}, 1:2, 1:2) with eltype Float64:
0.135318  0.958828
0.715183  0.136236

julia> sX = rand(2,2)
2×2 Matrix{Float64}:
0.57039   0.679983
0.914289  0.189875

julia> X
3×3 Matrix{Float64}:
0.135318  0.958828  0.0860895
0.715183  0.136236  0.542946
0.513489  0.974295  0.155925

As we can see, the attempt to assign to ‘sX’ creates a copy after all. Doesn’t this behavior annul the purpose of creating a view?

‘sX += I(2)’ would create a new array in similar fashion.

Only an additional slicing of ‘sX’ can change the values of the original array.

julia> sX[:] = rand(2,2)
2×2 Matrix{Float64}:
0.349409   0.523326
0.0127737  0.0894792

julia> X
3×3 Matrix{Float64}:
0.349409   0.523326   0.0860895
0.0127737  0.0894792  0.542946
0.513489   0.974295   0.155925
2

is an assignment to the variable sX, not the elements of the view previously assigned to it. Slicing specifies elements, so would sX .= rand(2,2).

5 Likes
3

= doesn’t mean what you think it means. See “assignment vs. mutation” in the Julia manual.

To replace the contents of an array with another array, you can do a[:] = … (setindex!) or a .= … (in-place broadcasting), or a loop where you set elements individually, or some in-place function call like rand!(a).

3 Likes
4

sXis a view of an array.
I don’t understand what is going on under the hood here.

sX += Ishould not produce a copy, because it has implicitly changed the type of sXfrom SubArray to Array.
I see only one explanations that would explain this behavior,

namely that there is promotion defined for SubArrayand Array, which is being called here

but I would argue that it should not behave like this:

5

The confusion has nothing to do with views or even with arrays. It’s just that, as indicated above, you misunderstood the semantics of the = syntax:

  • a = b: assign the value b to variable a

  • a[i] = b: setindex!(a, b, i)

  • a += b: a = a + b

1 Like
6

Forget what sX was.
Whatever the name sX was referring to, with sX = something_else you are telling Julia that the name sX must now point to the new something_else object, it doesn’t change in any way the original object sX was referring to.

2 Likes
7

I’m starting to suspect DanimirD is expecting behavior found in another language, which isn’t a good assumption to make unless specified otherwise e.g. Octave from MATLAB. In any case, it’ll help to speak more precisely than this to avoid confusion from other backgrounds.

The very first sentence of the Manual page Variables states plainly:

A variable, in Julia, is a name associated (or bound) to a value.

A variable isn’t a named value as you might find in other languages, it’s really just a name. This is a common paradigm in dynamically typed languages; if you allow variables to have any type at runtime, then it doesn’t make sense to make a value or its type part of the variable’s identity on the language level (though compiler implementations can leverage some patterns). So a = b really means:

  1. access the value assigned to variable b
  2. assign said value to variable a

Likewise, sX = rand(2,2) means:

  1. access the function assigned to rand
  2. call said function for inputs (2,2) to get a 2x2 Matrix{Float64}
  3. assign said matrix to variable sX

Note that the view previously assigned to sX was not involved at all, in fact it just lost its only binding so it doesn’t exist anymore on the language level (runtime implementation would have to free that memory automatically).

We may say “value b”, “function rand”, or “view sX” as shorthand, especially when the variable is const in a global scope, but it always means “the value currently assigned to this variable.”

1 Like
8

I guess Danimir was wondering the last command in the following block

julia> X = rand(3,3);

julia> sX = view(X, 1:2, 1:2)
2×2 view(::Matrix{Float64}, 1:2, 1:2) with eltype Float64:
 0.0748389  0.712181
 0.247739   0.453418

julia> sX = sX + [1 0; 0 1]
2×2 Matrix{Float64}:
 1.07484   0.712181
 0.247739  1.45342

(Because I also do not fully understand)
What are the precise processes happened when julia perform
sX + [1 0; 0 1]?
the former is a view, the latter is a Matrix.
Is there any middleman object? or allocations? (I don’t fully understand the meaning of “allocation” in computer systems)

Perhaps this makes some sense

julia> const a = [1];

julia> const b = [1];

julia> @allocated a + b
64

julia> @allocated broadcast!(+, a, a, b)
0

julia> println("a = $a, b = $b")
a = [2], b = [1]

julia> @allocated broadcast!(+, a, a, b)
0

julia> println("a = $a, b = $b")
a = [3], b = [1]
9

The above code is equivalent to this

new_array = sX + [1 0; 0 1]
sX = new_array

So you create and allocate a new array, and then you assign this new array to the variable sX. The old view that sX pointed to is gone.

If you want to update the value that sX refers to instead, you write

sX .+= [1 0; 0 1]
# or
sX .+= I

The example in the OP is best handled with

sX .= rand.()

The distinction between assignment (=) and mutation (.=) is crucial.

4 Likes
10

I have difficulty understanding why is this syntax.

I can explain why it is not sX .= rand(), which is because rand() only once.

But I don’t understand why can you put that dot behind rand and before an empty ().

I don’t think this reads rational.

Is there a functional API for this that attains the same or better performance?

11

It works because multiple dots in an expression fuse into a single broadcast call over all elements. For every element in sX there is a call to rand, which is directly input in the corresponding array position.

An alternative syntax is rand!(sX), but that is not as general, since sX .= f.() will automatically work for any function f (with a compatible return type), while f!(sX) would need to be implemented specially.

1 Like
12

Thanks, very useful! I’ve bookmarked this usage.

I still don’t understand why there isn’t a functional API counterpart.
I tried some methods which doesn’t produce good results in terms of performance

import Random.seed! as seed!
N = 9999

seed!(1);
a = -rand(N, N);
seed!(2);
@time @. a = rand(); # 0.103707 seconds

seed!(1);
a = -rand(N, N);
seed!(2);
@time broadcast!(_ -> rand(), a, a); # 0.160028 seconds (150.04 k allocations: 7.617 MiB, 42.89% compilation time)

seed!(1);
a = -rand(N, N);
seed!(2);
@time foreach(i -> setindex!(a, rand(), i), eachindex(a)); # 3.865129 seconds (199.97 M allocations: 2.980 GiB, 19.73% gc time, 0.30% compilation time)

is the current sX .= f.() style (or @. sX = f()) very flexible then? I doubt this is something still to be determined.

My doubt is because, sometimes, I need to pass an address to an API:
for example, the last arg view(X, j) here is an address

Gurobi.GRBgetdblattrelement(o, "X", 
        Gurobi.c_column(o, JuMP.index(x[j])),
        view(X, j))

in my this post Efficiently Retrieving Variable Values after Gurobi optimization with JuMP - #23 by WalterMadelim.

In this circumstance I can’t figure out a way to use the highly efficient dot syntax.

13

I am not able to understand the Jump syntax. Maybe you could give a more self-contained example of what you need?

I think some of your timing examples are down to suboptimal benchmarking. Take a look at BenchmarkTools.jl. I don’t think any of your code snippets would allocate if you benchmark correctly.

2 Likes
14

Thank you for the elaborate answer.

Yes, I do indeed expect a view to behave more like a reference.

And furthermore I expect it to be type stable.

This is Julia behavior btw:

julia> sX::Int = 1
1

julia> sX += 7//3
ERROR: InexactError: Int64(10//3)

This is kind of what I would consider “good” behavior: Throw an error if the assignment to a “reference” is not possible with the given expression - and I would be fine with that; much more than with the actual behavior of views.

I don’t think that implicit copies and or/conversions of the “view”-type (or SubArray) are the intended behavior of a programmer who wants to avoid copies in the first place.

Thanks for the explanation, but in that case my conclusion is that views are very limited in their utility.

15

They indeed do not:

julia> @btime @. $a = rand();
  119.171 ms (0 allocations: 0 bytes)

julia> @btime broadcast!(_ -> rand(), $a, $a);
  120.437 ms (0 allocations: 0 bytes)

julia> @btime foreach(i -> setindex!($a, rand(), i), eachindex($a));
  115.082 ms (0 allocations: 0 bytes)
2 Likes
16

I’m pretty sure none of these things have anything to do with the behavior of views. Views are not converted to a different type, and you can easily use them without causing allocations.

It seems that you are are probably not that familiar with dynamically typed languages.

2 Likes
17

in this case you get an error message purely because the type annotation ::Int, this code is still equivalent to:

sX::Int = 1

temp = sX + 7//3

sX = temp

Which would translate to this for your original problem:

julia> @views sX::SubArray = X[1:2,1:2];

julia> sX = rand(2,2)
ERROR: MethodError: Cannot `convert` an object of type
  Matrix{Float64} to an object of type
  SubArray

what you want is:

julia> sX .= rand(2,2)
2×2 view(::Matrix{Float64}, 1:2, 1:2) with eltype Float64:
 0.491047  0.707516
 0.649978  0.503986
1 Like
18

or

sX .= rand.()

or

using Random
rand!(sX)

which don’t allocate extra memory.

1 Like
19

I’m fairly convinced by now that you were and still are thinking of Julia variables as named values in statically typed languages.

  1. Julia is not statically typed, it’s strictly dynamically typed.
  2. Type stability is a practice to leverage the primary implementation’s type inference and compiler optimizations. Those are not language semantics, in fact the language can be implemented without type inference or a compiler (JuliaInterpreter does this to a large degree). Assignments definitely do not associate the value’s type with the variable; as I said before, this goes against the variable-object model.
  3. Loosely speaking, a view or a variable can be called “references”, but they are not referring to values the same way, and definitely not the same way as language-level pointers or references in other languages. The primary implementation of Julia has a core in C so there are memory addresses in some objects, but pointers don’t exist on the Julia language level. It’s not possible to exactly translate such semantics from another language to Julia.

Despite how this looks, this isn’t static typing nor is it making any part of the value part of the variable’s characteristics. In fact, you can declare a type for the variable without assigning/defining it, though you wouldn’t want to access it until it has been:

julia> global y::Int

julia> y
ERROR: UndefVarError: `y` not defined in `Main`
Suggestion: add an appropriate import or assignment. This global was declared but not assigned.

The declared type doesn’t have to exactly match the assigned value’s type either. Type declarations on variables add type checks to the assignment and try to convert the value to the type (on the language level, but the compiler can optimize this away). So, an abstract type for the variable is perfectly compatible with values of concrete subtypes, and the value’s type selects methods.

julia> z::Number = 1
1

julia> typeof(z)
Int64

julia> foo(::Int) = "hello";

julia> foo(::Number) = "goodbye";

julia> foo(z)
"hello"

Then don’t write the operations that the language specifies to do so. This isn’t a language where += has a separate “in-place” implementation from +, a += b literally means a = a + b. Sure, the compiler optimizes to in-place operations for some ubiquitous primitive types like i += 1 incrementing an integer, but that’s irrelevant on the language level. So, sX += I(2) just means assigning sX + I(2) to sX. Adding a SubArray and Diagonal requires promotion to a compatible output type in any language, and like the explanation in my previous comment, the view is also uninvolved in this reassignment to the output Matrix. If you want in-place addition to the view assigned to sX instead of reassigning sX, then you have to write sX[:,:] += I(2) or sX .+= I(2) to specify the view’s elements. There’s no way around that, that’s just how this language works.

20

I feel somewhat strange to add a $ at the left-hand-side of a =.
Despite this, there is still one difference:

In my original example, the three a outcomes are identical. Meaning that the three methods are from the same starting point and attains the same destination.

But with your @btime they don’t. Due to sample number? I don’t know.

Edit: Okay I’ve sorted it out :blush:

import Random.seed! as seed!
using BenchmarkTools
function get_starting_a()
    seed!(1);
    -rand(N, N);
end;

julia> N, S, E = 999, 98, 99;

julia> a = get_starting_a(); seed!(2);

julia> @btime(foreach(i -> setindex!($a, rand(), i), eachindex($a)); samples=S, evals = E);
  930.328 μs (0 allocations: 0 bytes)

julia> a = get_starting_a(); seed!(2);

julia> @btime(@. $a = rand(); samples=S, evals = E);
  951.779 μs (0 allocations: 0 bytes)

julia> a = get_starting_a(); seed!(2);

julia> @btime(broadcast!(_ -> rand(), $a, $a); samples=S, evals = E);
  900.982 μs (0 allocations: 0 bytes)

It’s indeed surprising that now @. is the slowest🙃