# Column-wise operations on matrices, allocating array views

**URL:** <https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331>\
**Category:** General Usage\
**Created:** [November 27, 2017, 1:15pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331 "2017-11-27T13:15:11Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 27, 2017, 1:15pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/1 "2017-11-27T13:15:12Z")

</div>

Hi,  
I find myself needing to do column-wise operations on matrices a lot (i.e. on A[:,j]). Or, spoken differently, I would want to use my matrix almost like a Vector{SVector}, i.e. A[j][:], which has identical data layout.

Ok, what do I need to do? Well, all the things one would like to do with vectors: arithmetic, swap, pass to functions.

Broadcast and array views have nice syntax but don’t work (they allocate). Vector{SVector} has nice syntax (I need to specialize on the number of rows though, and remember the changed order of dimensions), but crucially does not give by-reference-passing to functions.

So, I wanted to ask what other people are using. Of course I could write helper functions that try to behave like the broadcast/arrayview but keep the array-ref and the index separate; then I get no allocations, at the price of code-duplication (on functions that need to be capable of accepting both arrays and array\_views).

I could also go for writing explicit loops in all instances, but this kinda defies the point of using a high-level language (you wouldn’t even do this in C).

I also could go for a horrible “C-style” view (pointer), which does not protect the underlying matrix from gc. Sprinkling gc\_preserve might make this marginally better.

So: How are you people dealing with this problem? Any kind of julia-array that already solves this?

PS.  
Aggressive elimination of allocations probably won’t solve this. It will be solved once immutables containing reference-fields become bitstype, for the sake of code\_native (allocation/arg-passing/array-storage).

See also [https://github.com/JuliaStats/Distances.jl/issues/83](https://github.com/JuliaStats/Distances.jl/issues/83).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 27, 2017, 2:10pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/2 "2017-11-27T14:10:54Z")

</div>

Use `Vector{Vector}` or `Vector{MVector}`?

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [November 27, 2017, 2:54pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/3 "2017-11-27T14:54:38Z")

</div>

> I also could go for a horrible “C-style” view (pointer), which does not protect the underlying matrix from gc

Not recommended, but certainly possible. I used that approach in [NNLS.jl](https://github.com/rdeits/NNLS.jl/blob/7017cbbedfb40a4dac52c448545d217f685e4150/src/NNLS.jl#L211) where it was helpful in porting some existing Fortran code.

But in general I would agree with @stevengj: is there a compelling reason to have a `Matrix{T}` rather than a `Vector{Vector{T}}`?

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [November 27, 2017, 2:56pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/4 "2017-11-27T14:56:08Z")

</div>

> [@foobar\_lv2](#):
>
> crucially does not give by-reference-passing to functions

Also, what exactly do you mean by this? Everything in Julia has pass-by-reference (or, I suppose, pass-by-pointer) semantics, including SVectors. They’re just immutable, so they _might_ be copied if the compiler thinks that’s helpful, but whether a copy is made is irrelevant to the semantics.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 27, 2017, 3:11pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/5 "2017-11-27T15:11:18Z")

</div>

> [@rdeits](#):
>
> is there a compelling reason to have a Matrix{T} rather than a Vector{Vector{T}}?

Cache, indirection, storage overhead. Matrix{T} is equivalent to Vector{SVector}. If the inner vector is very large, then then Vector{Vector{T}} or array views are both fine (because I pay a constant price, amortized over a very large column). If the inner vector is very small, then SVector is fine. But it would be nice to use the same code for both.

Cases where you want to store a million 10-dim datapoints are terrible for Vector{Vector}.

> [@rdeits](#):
>
> Also, what exactly do you mean by this? Everything in Julia has pass-by-reference (or, I suppose, pass-by-pointer) semantics, including SVectors.

really?

As in: I have A=Vector{SVector}, and want to evaluate f(A[j]). Now the compiler either needs to make a copy or needs to: (1) defend against me changing A[j] (by overwriting with a new SVector) and (2) needs to keep alive A (because A[j] is allocated somewhere inside of A’s arraydata).

If the f is inline, then llvm should be able to avoid making this copy; else this should not be nice.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 27, 2017, 8:47pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/6 "2017-11-27T20:47:50Z")

</div>

It’s a little better on 0.7, because the compiler is smarter about eliding things and inlining “cheap” functions. For example,

```julia
julia> function sumcols!(dest, A::AbstractMatrix)
           _, indc = indices(A)
           @assert indices(dest) == (indc,)
           for i in indc
               @inbounds dest[i] = mysum(view(A, :, i))
           end
           dest
       end
sumcols! (generic function with 1 method)

julia> function mysum(v)
           s = 0.0
           @inbounds for x in v
               s += x
           end
           s
       end
mysum (generic function with 1 method)

julia> A = rand(1000,1000);

julia> dest = Vector{Float64}(1000);

# After warmup
julia> @time sumcols!(dest, A);
  0.001325 seconds (4 allocations: 160 bytes)

```

on 0.7 (but on 0.6 it has 1k allocations). However, if you replace `mysum` with `sum` then it allocates, because `sum` can’t inline.

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [November 27, 2017, 11:09pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/7 "2017-11-27T23:09:58Z")

</div>

> because sum can’t inline

Could you explain why that is? I dug around in the implementation of `sum` through `mapreduce`, and I don’t see any intentional `@noinline`, so I assume it’s something more subtle that prevents `sum()` from inlining?

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 27, 2017, 11:26pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/8 "2017-11-27T23:26:10Z")

</div>

Most likely just size. We don’t inline functions if we estimate their runtime cost to be significantly higher than that of a function call, and `sum` is more complex than `mysum` so probably gets penalized more heavily. More detail [here](https://docs.julialang.org/en/latest/devdocs/inference/#The-inlining-algorithm-(inline_worthy)-1).

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [November 27, 2017, 11:41pm UTC](https://discourse.julialang.org/t/column-wise-operations-on-matrices-allocating-array-views/7331/9 "2017-11-27T23:41:25Z")

</div>

I see. Thank you!