R: normalize (mean=0, variance=1) function that can be used with apply() on arbitrary axis

Question

I'd like to have a function that normalizes an array along a given axis. Basically, what I want to be able to write is:

apply(X, axis, normalize)

and it should normalize my arrays by rows, columns, ..., or along the n-th dimension.

This sounds like a very common problem, yet R seems to be unable to solve it without having to use quirks.

Consider the following matrix and function:

> m = matrix(1:4,nrow=2,dimnames=list(c('a','b'),c('x','y')))
> m
  x y
a 1 3
b 2 4

normalize=function(X) {
    X = X - mean(X)
    X = X/sd(X)
    return(X)
}

If I use apply() with axis numbers, this is what comes out:

> apply(m, 2, normalize)
           x          y
a -0.7071068 -0.7071068
b  0.7071068  0.7071068

Here everything's fine

> apply(m, 1, normalize)
           a          b
x -0.7071068 -0.7071068
y  0.7071068  0.7071068

Here, the matrix has been transposed (this is what I don't want). I want to retain the original dimensions for n-dimensional arrays.

Possible answers that are not working are:

• Use t(), apply it across rows and then transpose back: this does not work on multi-dimensional arrays
• Use function x of package y: likely the same as above
• Write an if clause for each dimension: the solution should be generally applicable, not only for some hardcoded cases

---

edit: with Roland's suggestions, I will use the following function instead of apply:

array_apply = function(X, along, FUN) {
    X = as.array(X)
    ndim = c(1:length(dim(X)))

    preserveAxes = ndim[ndim != along]
    orderAxes = c(along, preserveAxes)

    X = apply(X, preserveAxes, FUN)
    return(aperm(X, orderAxes))
}

which behaves the following way:

> m = matrix(1:4,nrow=2,dimnames=list(c('a','b'),c('x','y')))

> array_apply(m,1,normalize)
       x          y
a -0.7071068 -0.7071068
b  0.7071068  0.7071068

> array_apply(m,2,normalize)
       x         y
a -0.7071068 0.7071068
b -0.7071068 0.7071068

Answer 1

Let's use a better example to see what happens:

m = matrix(c(1,2,7,10,12,18),nrow=2)
#     [,1] [,2] [,3]
#[1,]    1    7   12
#[2,]    2   10   18
apply(m, 1, scale)
#            [,1] [,2]
#[1,] -1.02888681   -1
#[2,]  0.06052275    0
#[3,]  0.96836405    1

As you see scale is applied to the matrix rows and the resulting vectors combined as columns of the result matrix.

apply(m, 2, scale)
#           [,1]       [,2]       [,3]
#[1,] -0.7071068 -0.7071068 -0.7071068
#[2,]  0.7071068  0.7071068  0.7071068

Now scale is applied to the matrix columns, but the resulting vectors are again combined as columns of the result matrix. And indeed this is documented in help("apply"):

If each call to FUN returns a vector of length n, then apply returns an array of dimension c(n, dim(X)[MARGIN]) if n > 1.

Edit:

You can avoid the problem if you do something like this:

(m-apply(m,1,mean))/apply(m,1,sd)
          [,1]       [,2]      [,3]
[1,] -1.028887 0.06052275 0.9683641
[2,] -1.000000 0.00000000 1.0000000