Re: Vector approach to row margin frequencies

[Jaroslav Hajek]

On Thu, Jun 25, 2009 at 9:46 AM, Judd Storrs<address@hidden> wrote:
> On Thu, Jun 25, 2009 at 12:25 AM, Jaroslav Hajek <address@hidden> wrote:
>>
>> Surely is, because reshape is an O(1) operation (shares data, only
>> changes their interpretation). However, reshape/diag/reshape is only
>> sufficient if you want to scale the leading/trailing dimension of a
>> N-d array. If you need to scale a dimension in the middle (which is
>> uncommon, but can happen), you can use permute it to move it to the
>> end (or beginning, but end is better), or just turn again to repmat or
>> bsxfun.
>
> Don't underestimate diag() for the middle dimension ;-)
>
> X = ones(256,256,256) ;
> f = hamming(256) ;
>
> # merge dimensions 2 and 3, use diag
> tic
> g = diag(repmat(f,256,1)) ;
> reshape(reshape(X,256,65536)*g,256,256,256) ;
> toc
> Elapsed time is 0.150348 seconds.
>
> # merge dimensions 1 and 2, use diag
> tic
> g = diag(reshape(repmat(f',1,256),65536,1)) ;
> reshape(g*reshape(X,65536,256),256,256,256) ;
> toc
> Elapsed time is 0.1631 seconds.

Aww, cute. Manipulating both g and X to avoid permuting dimensions of
X didn't occur to me. This is of course much better than my permute()
solution. FYI, the permute() solution it seems 3x slower than the
above.

> # repmat and perform 3D element-by-element multiplication
> # most of the performance loss is in creating the full 3D factor
> tic
> g = repmat(f',[256,1,256]) ;
> X.*g ;
> toc
> Elapsed time is 0.318903 seconds.
>
> # bsxfun
> # Performance loss is presumably due to repeated function evaluations
> tic
> g = reshape(f,1,256,1) ;
> bsxfun(@times,X,g) ;
> toc
> Elapsed time is 2.264 seconds.

I'd thought bsxfun would be the loser. I intend to put in some
optimizations in the future.

> Is there any reason why bsxfun style singleton dimension expansion shouldn't
> be the *default* behavior for the builtin element-wise operators (.*, .+,
> .-) in octave? bsxfun also seems like the logical extension of the scalar
> with matrix cases.
>
> Matlab compatibility could be a reason, but I don't think it would break any
> working Matlab code (unless perhaps some strange code is relying on the
> error handling).  I think the gain in expressivity and simplicity for octave
> code could very well be worth it?
>
> Getting back to the original problem about scaling the matrix by row or by
> column sums if "./" behaved like bsxfun but avoided the repeated function
> evaluations then the blatently obvious expression of the algorithm would
> both work and be at least as fast as using diag,

Probably almost exactly as fast as diag, at least without repmat.

> while also scaling to
> higher dimensions sanely:
>
> X ./ sum(X,1)
> X ./ sum(X,2)
> X ./ sum(X,3)
> etc
>

I think the idea has been brought up multiple times in the past. A
major problem is implementation: somebody has to do it. Other minor
problems may include incorrect code that suddenly does some "magic"
rather than fail, and incompatibility with Matlab - some Octave users
limit themselves to Matlab-compatible subset. Scaling using diag, for
instance, has the advantage that it does work in Matlab.
But on general, I have no objections to that idea. However, in my
working, scaling/unscaling matrices is by far the most common
operation of this kind, and that is now solved elegantly using
diagonal matrices (which also work with sparse matrices efficiently),
so I have little motivation to actually do it.
I thought about optimizing bsxfun several times, but except the matrix
scaling case, It seems to me the result of bsxfun is usually a very
temporary expression that needs to be further reduced. A typical
instance is calculating the pairwise distances of two sets vectors. A
bsxfun + norm implementation would be better (with optimized bsxfun)
than a repmat + repmat + norm one, but still significantly slower than
a compiled function using loops.

-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz