What are the reasons for this benchmark result?

Question

Two functions that convert a rgb image to a gray scale image:

function rgb2gray_loop{T<:FloatingPoint}(A::Array{T,3})
  r,c = size(A)
  gray = similar(A,r,c)
  for i = 1:r
    for j = 1:c
      @inbounds gray[i,j] = 0.299*A[i,j,1] + 0.587*A[i,j,2] + 0.114 *A[i,j,3]
    end
  end
  return gray
end

And:

function rgb2gray_vec{T<:FloatingPoint}(A::Array{T,3})
  gray = similar(A,size(A)[1:2]...)
  gray = 0.299*A[:,:,1] + 0.587*A[:,:,2] + 0.114 *A[:,:,3]
  return gray
end

The first one is using loops, while the second one uses vectorization.

When benchmarking them (with the Benchmark package) I get the following results for different sized input images (f1 is the loop version, f2 the vectorized version):

A = rand(50,50,3):

| Row | Function | Average     | Relative | Replications |
|-----|----------|-------------|----------|--------------|
| 1   | "f1"     | 3.23746e-5  | 1.0      | 1000         |
| 2   | "f2"     | 0.000160214 | 4.94875  | 1000         |

A = rand(500,500,3):

| Row | Function | Average    | Relative | Replications |
|-----|----------|------------|----------|--------------|
| 1   | "f1"     | 0.00783007 | 1.0      | 100          |
| 2   | "f2"     | 0.0153099  | 1.95527  | 100          |

A = rand(5000,5000,3):

| Row | Function | Average  | Relative | Replications |
|-----|----------|----------|----------|--------------|
| 1   | "f1"     | 1.60534  | 2.56553  | 10           |
| 2   | "f2"     | 0.625734 | 1.0      | 10           |

I expected one function to be faster than the other (maybe f1 because of the inbounds macro).

But I can't explain, why the vectorized version gets faster for larger images. Why is that?

Answer 1

The answer for the results is that multidimensional arrays in Julia are stored in column-major order. See Julias Memory Order.

Fixed looped version, regarding column-major-order (inner and outer loop variables swapped):

function rgb2gray_loop{T<:FloatingPoint}(A::Array{T,3})
  r,c = size(A)
  gray = similar(A,r,c)
  for j = 1:c
    for i = 1:r
      @inbounds gray[i,j] = 0.299*A[i,j,1] + 0.587*A[i,j,2] + 0.114 *A[i,j,3]
    end
  end
  return gray
end

New results for A = rand(5000,5000,3):

| Row | Function | Average  | Relative | Replications |
|-----|----------|----------|----------|--------------|
| 1   | "f1"     | 0.107275 | 1.0      | 10           |
| 2   | "f2"     | 0.646872 | 6.03004  | 10           |

And the results for smaller Arrays:

A = rand(500,500,3):

| Row | Function | Average    | Relative | Replications |
|-----|----------|------------|----------|--------------|
| 1   | "f1"     | 0.00236405 | 1.0      | 100          |
| 2   | "f2"     | 0.0207249  | 8.76671  | 100          |

A = rand(50,50,3):

| Row | Function | Average     | Relative | Replications |
|-----|----------|-------------|----------|--------------|
| 1   | "f1"     | 4.29321e-5  | 1.0      | 1000         |
| 2   | "f2"     | 0.000224518 | 5.22961  | 1000         |