mean(X; dims=2) not the same as mean of each row

[sebastianskejoe]

Not sure if this is expected and/or acceptable behaviour, but I feel that it should at least be documented that mean(X; dims=2) can differ from [mean(X[i,:]) for i in 1:size(X,1)] - at least if the means are very small. I believe this is due to mean(X; dims=n) using mean! to calculate the mean over a given dimension, where as taking the mean of row calculates a "regular" mean (sum of elements divided by number of elements). Example output:

julia> [mean(X[i,:]) for i in 1:size(X,1)]
12-element Array{Float64,1}:
 -8.116963891148366e-16
 -1.64991477270683e-16
  1.450999814128156e-16
 -1.8621292779346593e-16
 -3.096288657565714e-16
 -9.86864910777917e-18
 -9.945747928933694e-18
 -9.398346298736568e-17
 -1.5034270125132327e-17
  6.908054375445419e-17
 -1.7578531223231644e-17
 -4.271274691960672e-17

julia> mean(X; dims=2)
12×1 Array{Float64,2}:
 -8.388351741612295e-16
 -1.8472877548624133e-16
  1.5743579279753957e-16
 -1.973151580397175e-16
 -3.367676508029642e-16
 -2.4671622769447924e-18
 -1.117932906740609e-17
 -9.274988184889329e-17
 -1.5651060694368526e-17
  6.661338147750939e-17
 -2.3746436915593628e-17
 -4.394632805807912e-17

[nalimilan]

Yes AFAIK that's because simple mean uses pairwise summation while mean over dimensions uses a plain mean for performance reasons. The docs should probably mention this.

[StefanKarpinski]

It seems like we could just use the same algorithm for both, no?

[nalimilan]

Not sure. When accumulating over non-contiguous slices, we need to store the result in an array (computing the sum one slice at a time would be much slower). So if we wanted to compute pairwise sum, we would have to allocate one such array per block.

[StefanKarpinski]

Note that there's no accumulation here, it's only reduction. The core of the sum reduction algorithm is:

https://github.com/JuliaLang/julia/blob/f0017a4964f47d515492167d/base/reduce.jl#L153-L176

I don't think any temporary arrays are required to make such a change.

[nalimilan]

Yes I meant reduction. But since reduction over dimensions computes the result for several slices in parallel (except for the special case of the first dimension), I think a temporary storage is needed for each block. See

julia/base/reducedim.jl

Lines 261 to 267 in f0017a4

	@inbounds for IA in CartesianIndices(indsAt)
	IR = Broadcast.newindex(IA, keep, Idefault)
	@simd for i in axes(A, 1)
	R[i,IR] = op(R[i,IR], f(A[i,IA]))
	end
	end
	end